Next Article in Journal
Analysis of the Expression Patterns of 13 DREB Family Genes Related to Cone-Setting Genes in Hybrid Larch (Larix kaempferi × Larix olgensis)
Next Article in Special Issue
Sliding Cutting and Cutting Parameters of Concentric Curvilineal Edge Sliding Cutter for Caragana korshinskii (C.K.) Branches
Previous Article in Journal
Water Use Efficiency of Five Tree Species and Its Relationships with Leaf Nutrients in a Subtropical Broad-Leaf Evergreen Forest of Southern China
Previous Article in Special Issue
Remote Monitoring of Amur Tigers in Forest Ecosystems Using Improved YOLOX Algorithm
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Deformable Shape Model for Automatic and Real-Time Dendrometry

Department of Forest Engineering, Resources and Management, College of Forestry, Oregon State University, Corvallis, OR 97331, USA
*
Author to whom correspondence should be addressed.
Current address: Silvx Labs, Missoula, MT 59802, USA.
Forests 2023, 14(12), 2299; https://doi.org/10.3390/f14122299
Submission received: 9 October 2023 / Revised: 11 November 2023 / Accepted: 21 November 2023 / Published: 23 November 2023
(This article belongs to the Special Issue New Development of Smart Forestry: Machine and Automation)

Abstract

:
We present a stereo image-based algorithm for tree stem diameter measurement and form analysis. The algorithm uses planar parametric curves to represent two-dimensional projections of tree stems in stereo images. The curves evolve according to an energy formulation based on the gradients of the images and inductive priors related to biomechanics and morphology of tree stems. After energy minimization, the curves are reconstructed to three dimensions, allowing for diameter measurements at any point along the height of the stem. We describe the algorithm and report the validation test results comparing predicted diameter measurements to external observations. Our findings demonstrate that the algorithm can automatically estimate diameters for trees within 20 m of the camera with an error of 5.52%. We highlight how this method can aid product value optimization through taper analysis and sweep or crook detection. A run-time analysis shows that the algorithm can estimate dendrometric variables for ten trees simultaneously at 15 frames per second on a consumer-grade computer. Furthermore, we discuss the opportunity to produce training data for machine learning algorithms that generalize across domains and eliminate the need to manually tune parameters.

1. Introduction

In forestry, decision-making often hinges on measurements obtained from individual trees. Given that such decisions can substantially influence both environmental and economic outcomes at broad scales, minimizing measurement errors and biases, and securing ample data to sufficiently decrease estimation errors are essential. Yet, due to the vast spatial reach of forests, the cost-effective collection of measurements remains a formidable challenge, making the reconciliation of these conflicting demands a central issue in forest inventory.
The emergence of contemporary sensor technologies, including digital photography and laser ranging, has considerably augmented our capacity to rapidly amass large data sets, often at decreased costs. These technologies have been employed in aerial surveys aimed at quantifying forest structure [1], estimating attributes [2,3,4], and mapping stems [5,6]. While aerial surveys provide extensive spatial information, they often depend on allometric correlations between sub-canopy attributes, such as diameter, and observations like tree height and crown width. However, allometric models utilizing aerial observations may only yield accurate predictions for select sub-canopy dendrometric variables of interest, falling short in providing data related to under-story density and composition. Consequently, there is a burgeoning interest in harnessing sensor technologies in terrestrial surveys to augment data collection efficiency and consistency.
One of the earliest instances of terrestrial remote sensing applications in dendrometry was introduced by Shelbourne and Namkoong [7]. They utilized two photographs of a tree stem, captured at 90-degree angles, to compute stem straightness in three dimensions (3D) by using the center-line of the two-dimensional (2D) projections. This approach was later refined by Hapca et al. [8], who integrated digital photography and extended the method to include stem taper profiles. They showcased the potential to automatically classify 3D stem reconstructions into form classes using their method [9].
Presently, terrestrial remote sensing research is primarily focused on the application of light detection and ranging (LiDAR) systems. Researchers have established the potential for constructing highly precise and intricate tree models [10,11] and automating stem detection and measurement [12,13,14] using LiDAR. However, despite the efficacy and accuracy of LiDAR systems in various forestry applications, their high cost and data processing demands render them impractical for certain uses.
The recent development of sophisticated computer vision algorithms and the availability of publicly accessible software has incited renewed research interest in photogrammetric approaches to dendrometry, specifically in the application of multi-view stereo and structure-from-motion (MVS-SFM) algorithms. This method involves converting a sequence of 2D images into a 3D point cloud, akin to data captured by a LiDAR sensor. The initial demonstrations of MVS-SFM application to dendrometry by [15,16] utilized a camera mounted on an unmanned aerial system (UAV) to acquire images of tree stems from various perspectives, reconstructing 3D point clouds and surface models. Subsequent studies by [17,18], using handheld cameras, discovered that diameter at breast height (DBH) estimations from MVS-SFM point clouds matched those estimated using terrestrial LiDAR. Later investigations into MVS-SFM-derived point clouds centered on detailed accuracy assessments [19], amalgamating aerial and terrestrial data [20], multi-camera configurations [21], and taper modeling [22].
To date, all applications of MVS-SFM for dendrometry follow a general workflow that includes data acquisition, offline processing for point cloud generation, and analysis of the point cloud through manual or automated measurements. While this approach offers satisfactory accuracy for most forestry applications and requires relatively inexpensive equipment, the computational expense of offline processing for point cloud generation can be significant, often taking hours or even days to complete. As such, in scenarios where immediate stem measurements are required, MVS-SFM may prove impractical. One such instance is markless prescriptions, wherein a description of the silvicultural treatment is provided to the operator of a harvesting machine tasked with selecting the trees to be cut. Such prescriptions frequently specify diameter limits, species preferences, and spacing requirements. In this context, a terrestrial remote sensing system that can deliver real-time, reliable measurements of trees within close proximity could prove advantageous for operators.
This paper presents a novel method for photogrammetric dendrometry, differing significantly from recent multi-view stereo and structure-from-motion (MVS-SFM) approaches. Our algorithm, based on fundamental photogrammetric principles, provides 3D reconstructions of tree stems and enables quantification of diameter, taper, sweep, and crook. This approach echoes earlier methods for analyzing tree stems in photographs [7] and incorporates the core principles of the active contour model presented in [23]. The inputs to the algorithm are the images from a calibrated stereo camera and the output of a bounding box object detector, trained on tree stems [24]. The algorithm automatically localizes the occluding contours of the tree stems within the bounding box in both the left and right stereo images, reconstructs the stems, and maps their positions to the object coordinate frame. Running in real-time on a modest graphics processing unit (GPU), the algorithm can provide reliable measurements for trees up to 20 m from the camera.
Aside from its direct applications in dendrometry, the approach presented here can also facilitate the creation of a training dataset for machine learning (ML) algorithms capable of segmenting tree stems in digital images. Although our method effectively localizes occluding contours in stereo images, the parameters of the algorithm require manual tuning specific to the problem domain. An ML-based approach, such as a region-based convolutional neural network (RCNN) [25], could offer a more generalizable and scalable solution. By using the 2D segmentation masks generated by our algorithm as training data for an ML-based segmentation model, we can develop a more robust system capable of segmenting tree stems under various environmental conditions and among different tree species. Consequently, this paper does not only present a novel photogrammetric approach to dendrometry but also lays the groundwork for future research that can integrate our work with state-of-the-art ML algorithms.
In this paper, we provide a comprehensive description of the algorithm, detailing the model, energy formulation, and optimization technique. We also delve into various algorithmic properties, showcase the results of a diameter measurement validation study, and discuss potential approaches to assess stem form and optimize product value from wood property inference. Figure 1 visually illustrates the major components of our proposed method.

2. Algorithm Description

We are given a stereo-rectified and row-aligned image pair, { I L ( x , y ) , I R ( x , y ) } , acquired from a stereo camera in the fronto-parallel configuration. Images are represented in function notation, I : Ω R , where Ω R 2 is the set of pixel coordinates in the image plane. We are also given the output of an object detector trained to detect tree stems [24]. The object detector outputs n bounding boxes, { Φ 1 , Φ 2 , , Φ n } , where each Φ i Ω . For example, we use the notation I Φ i : Φ i R to specify an image function within the domain of some arbitrary bounding box and I : Ω R to specify the entire image. We use upright, bold Greek letters to represent vector-valued functions in R 2 , e.g., χ : R R 2 , bold capital letters for matrices, e.g., A , and lowercase bold letters for vectors, e.g., x . We use the notation · as shorthand for · 2 , the Euclidean norm. Other notational conventions are apparent in context.

2.1. Deformable Stem Model

We use a set of parametric curves to represent the occluding contours of a tree stem projected on the image plane. We assume that there exists a curve of symmetry such that for all points on the curve, the unit normal, multiplied by some scalar, i.e., the radius, can be rotated about the tangent axis to the edge of the tree. Such a curve does not exist in 3D space, as the cross-section of a tree stem can only be approximated by a circle. However, this assumption holds for the two-dimensional case, and we exploit this symmetry to construct a shape model for tree stems.
Consider two curves, { λ : R R 2 , ρ : R R 2 } , representing the left and right contours of a tree stem projected in the x , y plane. We can construct these contours using the curve of symmetry, χ : R R 2 , and a taper function, τ : R R + , representing the change in radius along the stem. The left and right contours are defined as
λ ( s ) = χ ( s ) τ ( s ) n ^ ( s ) ,
ρ ( s ) = χ ( s ) + τ ( s ) n ^ ( s ) .
where s is the parameterization, s [ 0 , 1 ] , and the unit tangent and unit normal vectors of χ ( s ) are given by
t ^ ( s ) = χ ( s ) χ ( s ) ,
n ^ ( s ) = t ^ ( s ) t ^ ( s ) ,
Modeling the curves as presented requires an arbitrary number of parameters that depend on the level of discretization. In practice, discretization is often performed at the pixel level, resulting in hundreds of parameters that significantly increase computational demands during iterative optimization presented in Section 2.5. To address this challenge, we introduce a layer of abstraction to model tree stems that allows for efficient and compact representations without sacrificing the generality of the approach, even for instances of highly varied forms.
Let X = { ( x i , y i ) } i = 1 n denote an ordered set of coordinates in the plane. We use vectors to describe first and second entries independently, i.e., x = x 1 , x 2 , , x n T and y = y 1 , y 2 , , y n T . The coordinates in X correspond to n points on the curve of symmetry describing a tree stem. We fix the elements of y to be equally spaced along the height of the stem by
y i = i 1 n 1 h i = { 1 , 2 , , n } ,
where h is the height of the stem, equivalent to the number of pixel rows in the bounding box. We treat the vector x as a free parameter that can take on arbitrary values to represent an instance of the model. A curve can be constructed by interpolating X with a cubic spline interpolant S ( t ) C 2 [ y 1 , y n ] of the form
S ( t ) = S 1 ( t ) , y 1 t y 2 S 2 ( t ) , y 2 t y 3 S n 1 ( t ) , y n 1 t y n ,
where each S i ( t ) is a polynomial of degree three. We require that the interpolant passes through x , i.e., S ( y i ) = x i for all { i } 1 n , and constrain S ( y 1 ) = S ( y n ) = 0 to specify natural boundary conditions. The coefficients of the polynomials can be solved in O ( n ) using the tridiagonal matrix algorithm [27]. Using S x ( t ) to denote the spline function that depends on the parameter vector x , the curve of symmetry can now be defined as
χ ( s ) = S x ( h s ) , h s T ,
where s is the parameterization, s [ 0 , 1 ] .
Now that the curve of symmetry is specified, we define the taper curve in a similar manner. Let R = { ( r i , y i ) } i = 1 n be a set of n coordinates where the first entry elements, r = r 1 , r 2 , , r n T , specifies the radius of the stem at each corresponding second entry y , and each y i is defined by Equation (3). Using R we construct a cubic spline, S r ( t ) , satisfying the same conditions as for S x ( t ) , but in terms of the free parameter r . Following from Equation (5), a taper function describing the change in radius along the length of the stem is given by
τ ( s ) = S r ( h s ) , s [ 0 , 1 ] .
Computing the unit tangent and unit normal vectors as shown in Equation (2) with respect to χ ( s ) , we can use Equation (1) to define the curves representing the left and right contours of the stem, λ ( s ) and ρ ( s ) ; however, now they depend only on the parameter set { x , r } .
Figure 2 illustrates five instances of the Deformable Stem Model. Each instance is constructed using a different number of control points, n. As the number of control points increases, the number of parameters required to specify the model increases linearly by 2 n . Although an instance with a larger n can represent more complex stem forms, all stems we encountered during experimentation were sufficiently represented with 5–7 control points.

Extension to 3D

To extend our model to the object coordinate frame, R 3 , we can construct two sets of curves in stereo such that one corresponds to the left camera and the second to the right. We use superscripts L and R on the functions to denote a curve in the left and right image, respectively. Assuming that the stereo pair is rectified and row-aligned, the curves can be back-projected to the object coordinate frame via triangulation. First, we compute the inverse disparity, δ : R R + , given by
δ ( s ) = | χ x L ( s ) χ x R ( s ) | 1 ,
where the subscript on the functions index the corresponding position in the output vector. Using the inverse disparity and the stereo calibration parameters, the curve of symmetry in the left camera coordinate frame is representing in object coordinates by the function X : R R 3 , defined as
X ( s ) = χ x L ( s ) c x δ ( s ) b χ y L ( s ) c y δ ( s ) b f δ ( s ) b ,
where ( c x , c y ) is the principal point of the left image, f is the focal length in pixels and b is the baseline distance between the optical axes of the left and right cameras. The taper function is scaled to object coordinates by
T ( s ) = τ L ( s ) δ ( s ) b ,
where T : R R + . Given the back-projected curve of symmetry, X ( s ) , and the scaled taper function, T ( s ) , we can define a parametric surface function, Λ : R 2 R 3 , based on Rodrigues’ rotation formula to reconstruct the tree stem to 3D:
Λ ( s , θ ) = X ( s ) + T ( s ) I + sin θ t ^ ( s ) × + 1 cos θ t ^ ( s ) × 2 n ^ ( s ) ,
where s and θ are the parameterizations, s [ 0 , 1 ] and θ [ 0 , 2 π ] , I is a 3 × 3 identity matrix, t ^ : R R 3 and n ^ : R R 3 are the unit tangent and normal vector functions of X ( s ) , and [ t ^ ( s ) ] × is a skew symmetric matrix of the unit tangent vector defined as
t ^ ( s ) × = 0 t ^ z ( s ) t ^ y ( s ) t ^ z ( s ) 0 t ^ x ( s ) t ^ y ( s ) t ^ x ( s ) 0 .

2.2. External Energy

Active contour models use an external energy to attract contours to desirable features in the image. External energy is derived from the input image in such a way that minimum values correspond to the desirable features. Our objective is to localize the edges of tree stems, so we will formulate the external energy to be minimal near vertical edges in the input image. We also introduce a secondary external energy to add robustness to local minima during curve evolution.

2.2.1. Edge Map and Gradient Flow

We define an edge map I x : Φ i R + by differentiating the image I : Φ i R w.r.t x and taking the absolute value,
I x ( x , y ) = | I x | .
In practice, it is recommended to denoise the image prior to computing gradients via Gaussian smoothing or bilateral filtering. We avoid computing the gradient vector I (differentiating I w.r.t x and y) as edges that represent the boundaries of the stem will only have strong gradients in the x direction while gradients in the y direction could introduce spurious forces that adversely influence the convergence of the contours.
Differentiating the edge map, I x ( x , y ) , w.r.t x gives a function where the zero-crossings correspond to edges in the image. If we consider this as a functional taking a parametric curve as an input, then the sign and magnitude of the scalar values returned by the functional indicate the direction and speed in which to move the curve. This, however, assumes that the curve, at some point during evolution, is within the capture range of a boundary in the edge map, which is generally an unrealistic assumption. Xu and Prince [28] present an alternative formulation for external energy, called gradient vector flow (GVF), that overcomes issues with insufficient capture range by filling in regions where gradients are not present, effectively increasing the capture range of the edge map.
GVF is a vector field v ( x , y ) = u ( x , y ) , v ( x , y ) T that minimizes the energy
E u v = η ( u x 2 + u y 2 + v x 2 + v y 2 ) + I x 2 v I x 2 d x d y ,
where η is a regularization parameter that controls the relative importance of the first and second terms in the integrand. Since the Deformable Stem Model considers updates only in the x direction, we simplify the expression by omitting the y-component of the vector field yielding
E u = η ( u x 2 + u y 2 ) + I x x 2 u I x x 2 d x ,
where I x x is defined as
I x x ( x , y ) = I x x .
As shown by Xu and Prince [28], u can be found by expressing Equation (14) as a Euler equation,
η 2 u ( u I x x ) I x x 2 = 0 ,
where 2 is the Laplacian operator, and then solving via iteration by treating u as a function of time,
u t + 1 ( x , y ) = η 2 u t ( x , y ) u t ( x , y ) I x x ( x , y ) I x x ( x , y ) 2 .
Through experimentation, we found η = 0.2 to be a sufficient setting for the regularization parameter. We use the function G : R 2 R to represent the steady-state solution to u t ( x , y ) . We refer to G ( x , y ) as the gradient flow force function to indicate that it is not a vector field as originally proposed by Xu and Prince [28]. G ( x , y ) returns a signed scalar value where the sign denotes the direction of the force and the value denotes the magnitude of the force queried at the position variables x and y.

2.2.2. Contraction Energy

The edge map I x has the property that large values correspond to abrupt intensity changes in the horizontal direction. Ideally, large values would represent the boundaries of the stem while the rest of the domain exhibits small values. However, I x is likely to contain large values elsewhere; large values might exist in the background of the image caused by adjacent stems, branches, or other features in the image. We formulate a contraction energy based on the stereo disparity map to make the contours robust to local minima that do not correspond to gradients along the tree stem boundary. We compute an integer stereo disparity map for the entire image domain, D : Ω N , using semi-global matching [29]. We use a real-time GPU implementation of the semi-global matching algorithm presented by [30].
Let B : R 2 { 0 , 1 } be a logical function defined as
B ( x , y ) = 1 if d θ D Φ ( x , y ) d + θ 0 otherwise ,
where D Φ : Φ i N is the stereo disparity map within the domain of the bounding box, Φ i . The threshold value d can be readily found by computing the histogram of D Φ ( x , y ) and taking d as the argument maximum of the histogram. The range value θ is the expected range above and below d used to capture pixels of the tree stem projected on disparity planes other than d, typically 3 to 5. Next, we compute a row-wise distance transform, W : R 2 R , given by
W ( x , y ) = max u : B ( u , y ) = 1 1 2 ( u x ) 2 .
The form of Equation (19) becomes apparent when taking its derivative w.r.t. x,
W x ( x , y ) = W x .
W x ( x , y ) gives the signed distance in pixels along the x-axis to the closest pixel in B ( x , y ) = 1 . This is a desirable property since it guides the contours towards a coarse segmentation of the tree stem given by the disparity map making the evolution of the contours robust to local minima in the image.
We now have the necessary components to define the external energy functionals. For the left image we combine the negative of the edge map, I x ( x , y ) , so that strong edges correspond to minima, and the negative of the contraction energy, W ( x , y ) , so that the pixels where B ( x , y ) = 1 are minimum. We also scale the contraction energy by the scalar κ to balance it with the edge map so that the contraction energy becomes small as the contours approach the boundaries of B ( x , y ) . For the right image, we simply use the negative of the edge map since the contraction force is derived from the disparity map that corresponds to the left image frame. The stem model corresponding to the right image will be influenced by the contraction energy through the stereo constraints introduced in Section 2.4. Formally, the external energy for the left and right images are defined as
E ext L ( x , y ) = I x L ( x , y ) κ W ( x , y ) ,
E ext R ( x , y ) = I x R ( x , y ) .
The energy maps can be converted to horizontal force maps by taking their derivatives w.r.t. x. The gradient flow of the edge map, G ( x , y ) , defined in the previous section is used as the force map for edges while W x ( x , y ) is used for the contraction force,
F L ( x , y ) = G L ( x , y ) + κ W x ( x , y ) ,
F R ( x , y ) = G R ( x , y ) .
We found κ [ 0.05 , 0.1 ] to be good values for scaling the contraction energy.
Figure 3 illustrates the progression from an input image to the force maps. Note that the capture range in G ( x , y ) , Figure 3d, is far greater than the capture range in the force map I x x ( x , y ) , Figure 3c.

2.2.3. External Forces on the Deformable Stem Model

In this section we show how the external force map is used to calculate the external forces acting on the Deformable Stem Model. Recall the parameters x and r used to specify an instance of the model. These parameters are used to construct the left and right contours of a tree stem, λ ( s ) and ρ ( s ) , respectively. Since our external energy function, E ext ( x , y ) , is intended to be minimal at the boundary of the tree stem, our objective is to evolve the contours in such a way that each subsequent state results in a lower energy when integrated along the length of the contours. We can achieve this by subdividing the contours and integrating over the subintervals such that each subinterval corresponds to the control points in x and r . The integrals along the subintervals are averaged and placed in the force vectors f λ = f λ , 1 , f λ , 2 , , f λ , n T and f ρ = f ρ , 1 , f ρ , 2 , , f ρ , n T , then used to update x and r . This is similar to the approach taken by Wang and Boyer [31] in the Active Geometric Shape Model.
We use a = a 1 , a 2 , , a n + 1 T to represent the subintervals where each element is given by
a 1 = 0 ,
a i = i n 1 1 2 ( n 1 ) , { i } 2 n ,
a n + 1 = 1 .
The individual force components can be computed by integrating over the gradient flow force function along the sub-intervals for the left and right contours:
f λ , i ( j ) = 1 a i + 1 a i a i a i + 1 F j λ j ( s ) d s , { i } 1 n , j = { L , R } ,
f ρ , i ( j ) = 1 a i + 1 a i a i a i + 1 F j ρ j ( s ) d s , { i } 1 n , j = { L , R } ,
where, again, s [ 0 , 1 ] is the parameterization. The components of the force vectors can be interpreted as the average force acting on the subinterval of the curve. So, if the part of the curve bounded by the subinterval is positioned to the left of a boundary in the energy function, then the force will be greater than 0; a force less than 0 indicates that the part of the curve is positioned to the right of the boundary. A force of zero means that the curve segment is balanced on the boundary.
The parameter vector x , representing the control points for the curve of symmetry, can be updated to a state of lower energy by adding x and the vector sum of the forces acting on the left and right curves,
x t + 1 ( j ) = x t ( j ) + f λ ( j ) + f ρ ( j ) , j = { L , R } .
The parameter vector r , representing the control points for the taper curve, can be updated by adding the vector with the forces acting on the right contour less the forces acting on the left,
r t + 1 ( j ) = r t ( j ) + f ρ ( j ) f λ ( j ) , j = { L , R } .
To better understand these update equations, it is recommended that the reader uses Figure 3d,e to visualize different configurations of the left and right curves and how the update equations apply. Note that darker colors in the figures correspond to negative values while lighter colors correspond to positive values.

2.3. Internal Energy

We have yet to constrain the parameter set used to construct an instance of the Deformable Stem Model; the parameters x and r are controlled entirely by the external forces, some of which might not correctly represent the boundary of the tree stem. In this section, we introduce two biomechanically-inspired internal energy formulations that are simple analogies to how trees respond to physical forces in their environment.

2.3.1. Straightness Force

Trees have a natural tendency to grow straight. This is influenced by an internal control process called gravitropism, first recorded by Charles Darwin in 1880 [32]. Trees that grow against the gravity vector have the advantage that their weight is centered over the stem, thus the downward forces acting on the stem can be distributed. Trees also grow toward light sources through a weaker process, called phototropism, which can cause trees to lean or sweep toward openings in the canopy. Our objective here is to formulate an internal force for the stem model to encourage it to be straight while allowing for slight deviations to account for irregular stem forms.
We derive straightness energy using the normal equations for linear least squares regression. The energy of straightness corresponds to the sum of squared residuals from the best-fit line relating the parameter vector x and the arithmetic sequence y defined by Equation (3). Since the linear least squares model is undefined for a vertical line, we exchange the axes and define the design matrix on y by
Y = 1 y ,
and take x as the vector of dependent variables. The notation 1 denotes a vector of ones with the same dimension as y . Rearranging the normal equations for linear regression, we can get the predicted values for the best-fit line by
x ^ = Ξ x ,
where Ξ is the Hat matrix defined as
Ξ = Y ( Y T Y ) 1 Y T .
The Hat matrix can also be used to compute the residuals, ϵ , by subtracting it from the identity matrix and right multiplying with the dependent vector,
ϵ = I Ξ S x .
We take the Euclidean norm of the residual vector as the straightness energy,
E straightness = S x .
If the straightness energy is zero, then it is implied that the coordinates X , defined in Section 2.1, are co-linear.
In order for the straightness energy to interact with other forces in the model, it is necessary to consider the force components that make up the scalar quantity describing the energy. This can be accomplished by discretization of the individual forces in the time domain and taking small steps in the direction of the forces. We can subtract α S x from x , where α [ 0 , 1 ] , to transition x to a state with lower energy. For example, if α = 0 , then x α S x = x , and if α = 1 , then x α S x = Ξ x = x ^ . Using x t to index x at time t, we have the following update rule for minimizing the straightness energy in the time domain,
x t + 1 = x t α S x t .
Figure 4 illustrates an example of the internal straightness force acting on an instance of the Deformable Stem Model. The instance (a) has a straightness energy of 8.82. After 10 iterations with α = 0.1 the energy of the stem model shown in (c) is decreased to 3.07.

2.3.2. Monotonic Radii Force

In 1893, Metzger [33] provided an biomechanical explanation for the taper of tree stems. Metzger hypothesized that trees distribute radial growth along their stem to uniformly counterbalance the mechanical stress acting on the tree caused by lateral forces, e.g., wind. A beam of uniform resistance to these forces is a cubic paraboloid, a geometric shape that closely resembles that of a tree stem. In reality, trees have more complicated taper functions, as indicated in numerous studies [34,35,36]. Here, we enforce a simple constraint on the stem model to encourage the radius parameters to be monotonically decreasing with height to allow for variations among species and individuals.
Let A be a backwards difference matrix. We can construct this matrix by subtracting an n × n identity matrix, I n , by a super-diagonal shift matrix, U n , of the same shape,
A = U n I n .
If we take the matrix-vector product A r , then the positive values in the resulting vector indicate that the corresponding element in the radius vector, r , has a larger subsequent element, i.e., the projected width of the stem is increasing with height, which is generally not how trees taper. Negative values indicate that the width is decreasing, while zero indicates a constant taper. Our objective here is to formulate an energy that is zero when the radius vector of the stem is monotonically decreasing with height, and positive otherwise. We can achieve this by first passing the matrix-vector product A r through a linear rectifier to truncate all negative values to zero. The rectifier is a vector-valued function, Ψ : R n R n , defined as
Ψ ( p ) = max ( 0 , p 1 ) max ( 0 , p 2 ) max ( 0 , p n ) T ,
where p is a placeholder for the product of A r . Next, we average bidirectional differences of increasing values with the operation 1 2 A T Ψ ( A r ) . The sign and magnitude of the elements in the resulting vector guide the movement of the control points for the taper curve to a state of lower energy, where the energy, E taper , is given by
E taper = 1 2 A T Ψ ( A r ) ,
We can discretize the movement of r in the time domain by introducing the step size parameter, β , leading to the following update rule for minimizing taper energy:
r t + 1 = r t β 1 2 A T Ψ ( A r t ) .
Figure 5 shows the internal taper force acting on an instance of the Deformable Stem Model. The instance shown in (a) has an energy of 2.34. The resulting model following 10 iterations with β = 0.1 , instance (c) in the figure, has an energy of 0.72.

2.4. Stereo Energy

In this section we introduce the two final energy formulations. The energy terms presented here create a connection between the models in the left and right images. Specifically, information regarding the contours in one image is used to localize the contours in the other. This is especially useful when there are desirable features missing in one image that are present in the other.
We define a stereo smoothness constraint on the curves of symmetry in the left and right image as
E stereo-sym = 1 2 d χ L d x d χ R d x 2 .
This ensures that the disparities between the curves vary smoothly. E stereo-sym = 0 means that the derivatives of the curves w.r.t. x are exactly the same, while E stereo-sym > 0 means the disparities diverge somewhere along the length of the curves. Since χ ( s ) is a spline interpolation of the parameter vector x , we can define force update equations for the left and right model parameters as
x t + 1 ( L ) = x t ( L ) + μ D x t ( R ) D x t ( L ) ,
x t + 1 ( R ) = x t ( R ) + μ D x t ( L ) D x t ( R ) ,
where μ is the step size parameter and D is a finite difference matrix. We use the step size parameter to make a soft constraint on the differences. This allows the disparity between the curves to vary slightly when the tree stem is not exactly vertically aligned with the image plane.
The second stereo constraint forces the taper curve of the left and right model to be proportional,
E stereo-tap = τ L ( s ) γ τ R ( s )
where γ is some scalar that, when multiplied by one of the taper curves, minimizes their differences. Since the two taper profiles are projections of essentially the same face of the stem offset by a small baseline, we expect their taper functions to be identical. However, depending on the angle at which the stem is positioned w.r.t. the optical axes of the camera, the distance from the stem to the left camera can be different than the distance from the stem to the right camera. Thus, we allow the scale of the taper function to differ. This constraint can be converted to force update equations and applied to the left and right radius parameter vectors by first taking the average of the two vectors and then adding the sum of their difference divided by 2 n ,
r t + 1 ( L ) = 1 2 r t ( L ) + r t ( R ) + 1 2 n i n r t , i ( L ) r t , i ( R ) ,
r t + 1 ( R ) = 1 2 r t ( L ) + r t ( R ) + 1 2 n i n r t , i ( R ) r t , i ( L ) .

2.5. Optimization

In this section we show how we minimize the external, internal and stereo energy of the Deformable Stem Model via gradient descent. First, the parameters for the model, { x , r } , are initialized according to the width of the image domain Φ i . We set all values of x and r , for both the left and right models, equal to w 2 , where w is the number of pixel columns in Φ i . Thus, the left and right contours, when computed using Equation (1), are initialized at the left and right edges of the bounding box. We also initialize y with Equation (3) where h in the equation is set to the number of pixel rows in Φ i . The matrices S , D , and A are constants, so they are initialized and stored. The external force maps for the left and right images are computed using Equation (22).
Next, we concatenate all the force update equations for the left and right parameter vectors x ( L ) and x ( R ) ,
x t + 1 ( L ) = x t ( L ) α S x t ( L ) + μ D x t ( R ) D x t ( L ) + f γ ( L ) + f ρ ( L ) ,
x t + 1 ( R ) = x t ( R ) α S x t ( R ) internal + μ D x t ( L ) D x t ( R ) stereo + f γ ( R ) + f ρ ( R ) external ,
and for the parameter vectors r ( L ) and r ( R ) ,
r t + 1 ( L ) = r t ( L ) β 1 2 A T Ψ A r t ( L ) internal + 1 2 r t ( L ) + r t ( R ) + 1 2 n i n r t , i ( L ) r t , i ( R ) stereo + f ρ ( L ) f γ ( L ) external ,
r t + 1 ( R ) = r t ( R ) β 1 2 A T Ψ A r t ( R ) internal + 1 2 r t ( L ) + r t ( R ) + 1 2 n i n r t , i ( R ) r t , i ( L ) stereo + f ρ ( R ) f γ ( R ) external .
We iterate these equations and during each iteration we construct a cubic spline interpolant through { x ( L ) , r ( L ) } and { x ( R ) , r ( R ) } in order to compute the edge contours { λ L ( s ) , ρ L ( s ) } , and { λ R ( s ) , ρ R ( s ) } . This is necessary for calculating the external force vectors—the last terms in Equations (41) and (42). After each iteration, we check for convergence, where the convergence value is given by
ϵ = x t + 1 ( L ) x t ( L ) + r t + 1 ( L ) r t ( L ) + x t + 1 ( R ) x t ( R ) + r t + 1 ( R ) r t ( R ) .
The iteration is terminated when ϵ is less than some threshold, e.g., 0.1 . We do this to avoid issues with numerical instability when the parameters oscillate around a minimum. Once the model parameters have converged, we use the stereo calibration parameters to reconstruct the model using Equations (7)–(10).
Refer to Figure 6 for a diagram showing the flow of operations in the algorithm. We have referenced the core equations used in the text boxes between the visual representations. Figure 7 illustrates curve evolution during optimization of the Deformable Stem Model on six example trees. For each example, we show 5 snap-shots of the contours corresponding to 1, 5, 15, 30, and 50 iterations, respectively. The left-most image in each subfigure shows the position of the left and right edge contour following the first update after initialization based on the given bounding box. The right-most image shows the edge contours after convergence. The images in between are intermediate time steps. We also alternate left and right stereo images, so images 1, 3, 5 are from the left camera and 2 and 4 are from the right camera.

3. Application to Diameter Measurement

3.1. Materials

We tested the Deformable Stem Model using a custom-built stereo camera with a 32 cm baseline. The cameras were operated at VGA resolution ( 640 × 480 ) with a focal length of 3.0 mm. The cameras were calibrated for intrinsic and extrinsic parameters following methods presented in [37], and the images were stereo-rectified and row-aligned. We captured 7 stereo pairs within a sparse ponderosa pine (Pinus ponderosa Douglas ex Lawson) forest in western Montana, each photograph containing 2 to 4 trees within 20 m of the camera, for a total of 21 trees. We recorded a caliper measurement for each tree across the face of the stem pointing toward the camera at breast height (1.37 m) and the distance from the plane of the left camera imaging sensor to the center of each tree at breast height. Table 1 shows the mean, standard deviation, minimum, and maximum values for the collection data.
We also captured a second stereo pair at each plot after wrapping each tree with orange flagging at the height of measurement so that the exact point of measurement along the stem could be converted to a pixel location for validation purposes. We minimized the energy of the Deformable Stem Model on each tree using the parameter configuration shown in Table 2.

3.2. Analysis

First, we tested the accuracy of diameter measurements without using the camera calibration parameters to reconstruct the stem model. We do this to assess how potential errors in camera calibration contribute to the measurement accuracy. We use the notation τ ( y ) to denote the taper function at the observed breast height measurement position y. The taper function returns the number of pixels from the center of the tree to the edge, so scaling by 2 gives the number of pixels across the entire stem. The estimated diameter at breast height is given by
d ^ = 2 τ ( y ) z f ,
where z is the observed distance from the camera to the tree and f is the focal length. This implies that
d ^ 2 τ ( y ) z .
Thus, we can relate the observed diameter measurements, d, to the proportional estimates, 2 τ ( y ) z , without using the calibration parameter f. We use the linear regression parameters, β 0 and β 1 to compute the estimated diameter,
d ^ ( z , y ) = 2 τ ( y ) z β 0 β 1 .
Next, we estimate the diameter using the stereo disparity from the left and right stem models and the calibration parameters to estimate the distance to the tree stem. This estimate is not considered completely automatic since we use the observed height of measurement position y. The estimator for the distance is given by,
z ^ ( y ) = f b δ ( y ) ,
where b is the baseline distance between the cameras in meters and δ ( s ) is the inverse disparity function defined in Equation (7). Plugging z ^ ( y ) into Equation (44) yields
d ^ ( z ^ , y ) = 2 τ ( y ) z ^ ( y ) f .
Note that this equation can be simplified so that the variable f appears in the numerator and denominator. However, we keep the equation as is to indicate that z is a function of a specific height argument.
For the final diameter estimator, we replace the known height position variable, y, with an estimated measurement height. We used RANSAC [38] to extract the ground plane from the point cloud generated with the disparity map and used the best-fit plane as a reference to calculate breast height [39]. The estimator for the measurement height, y ^ , corresponding to the best-fit plane describing the ground, is defined as
y ^ = argmin s [ 0 , 1 ] | n ^ X ( s ) μ 1.37 | ,
where X ( s ) is the back-projected curve of symmetry defined in Equation (8), n ^ R 3 is the unit normal vector of the plane, μ R 3 is an arbitrary point on the plane and 1.37 is the distance above the ground to breast height in meters. The final diameter estimator is given by
d ^ ( z ^ , y ^ ) = 2 τ ( y ^ ) z ^ ( y ^ ) f .
For each estimator we calculated the root mean squared error (RMSE), defined as
1 N i N θ ^ θ 2 1 2 ,
the bias on the estimate, defined as
1 N i N θ ^ θ .
and the mean absolute error (MAE), given by
1 N i N | θ ^ θ | .

4. Results and Discussion

4.1. Diameter Measurement

Figure 8a illustrates the relationship between the ground truth diameter, d, and the projected width in pixels, 2 τ ( y ) , scaled by the ground truth distance, z. Figure 8b shows the relationship between z and z ^ ( · ) , where the open circles in the figure correspond to z ^ ( y ) . As the figure suggests, any existing calibration errors are minuscule. The estimated diameters according to Equation (50) are plotted against the observed diameters in Figure 8c. The residuals of the estimates are plotted against the observed distance in Figure 8d to show that the errors do not appear to be correlated with distance.
The RMSE, bias, and MAE of the estimates are presented in Table 3. We focus our discussion on RMSE as it penalizes large errors more than MAE. The RMSE for the calibration-free diameter estimate, d ^ ( z , y ) , was 1.61 cm with a small bias of 0.02 cm. We consider this estimate to have the highest potential accuracy among all other estimates assuming error-free distance observations. Using the camera calibration parameters to estimate the distance increased the RMSE of the diameter estimate by 0.15 cm and the bias by 0.11 cm. This marginal increase in error suggests that our camera calibration parameters are sufficiently accurate. The final diameter estimate, using both the estimated distance and measurement height, had an RMSE of 1.94 cm, 0.33 cm greater than the calibration-free estimator, and slightly overestimated the ground truth mean by 0.3 cm. In comparison to other recent studies that focus only on non-contact diameter measurement at breast height, our RMSE is within the range of their results: 7.9 cm [40], 2.71 cm [41], 1.28–2.57 cm [42], and 1.55 cm [43].

4.2. Run Time Analysis

The algorithm was implemented in the C++ programming language using the CUDA parallel computing platform for interfacing with parallelization elements of the GPU [44]. We performed a worst-case run time analysis on 25 images of trees ranging from 2–56 kilopixels (kp) in size. We considered this analysis a worst-case evaluation as we set the number of iterations for the optimization step to 100. This is approximately double the number of iterations required for convergence when we tested the algorithm on the 21 trees in the diameter validation experiment. In Table 4, we show the average time in milliseconds (ms) for the major components of the algorithm configured according to Table 2. The table also summarizes the time required for allocating GPU memory and uploading data to the GPU. These tests were performed using a NVIDIA GeForce GTX 780M GPU.
The average run time across all test images was 5.57 ms (179.53 frames per second). In the table, we do not include the run time of the stereo-matching algorithm that is required for computing the contraction force since it depends heavily on the choice of stereo-matching algorithm and implementation. The implementation we used [30], can perform stereo matching for a VGA resolution image in 11 ms (90 fps) on the GTX 780M GPU. This increases the run time of our algorithm to 16.57 ms and decreases the frame rate to 60 fps. Furthermore, this analysis summarizes the run time for a single tree. If we expect to have on average 10 tree detections at any given time, then the run time, including the time required for stereo matching, will increase to approximately 66 ms (15 fps).

5. Conclusions

In this paper, we demonstrated how 2D projections of tree stems in stereo images can be automatically reconstructed into 3D models in real time. We show how dendrometric attributes can be calculated from the models and provide an accuracy assessment of diameter measurements for 21 tree stems. The algorithm presented is based on the active contour model, a well-developed computer vision framework for boundary localization. The algorithm runs in real-time on a modest GPU and is intended to be used during forest harvesting operations for diameter measurement and stem form analysis of standing trees.
There are many attributes of a tree stem that influence its value. Stem size and shape are among the most important variables and can also be used to make optimal decisions regarding bucking positions and milling routines [45,46]. The algorithm presented in this paper can provide 3D stem models for product value optimization. Although our approach only provides caliper measurements for the sub-canopy portion of the stem, the shape of the stem is calculated in 3D. Since these models can be generated in real-time, they can be used in harvesting operations to make optimal bucking decisions and to calculate the value of the stem.
Research has also been conducted to relate standing tree shape to the spatial distribution of internal wood characteristics. Constant et al. [47] presented an analysis of poplar trees that suggests a strong relationship between local tree slope and the position of the pith in relation to the center of the stem. They also found that the concentration of reaction wood is related to the local slope of the stem. Since reaction wood and the eccentricity of annual rings are directly related to the quality and price of the products derived from trees, 3D models like those provided by our algorithm can be useful in determining the value of logs prior to harvesting.
The data used in our diameter validation experiment were collected in a park-like forest where the trees were uniformly distributed and similar in age class. We have also tested our algorithm in a natural forest where trees are clustered and vary in age class. Although our algorithm performs well in natural forests, it is unclear what density of trees or level of understory vegetation will limit the performance of the algorithm.
A limitation of this work is that it is based on only two unique views of the stem as opposed to many in the case of MVS-SFM, essentially providing caliper measurements of the stem. We do not, however, consider this to be a major limitation as caliper and d-tape measurements have been shown to be statistically similar [48]. Another limitation is that our method does not provide tree height, which is an important dendrometric variable. We recommend that future research in this arena consider the use of fish-eye or hemispherical lenses to capture the upper portions of tree stems.
The analysis presented in this paper did not provide insight to the maximum distance a tree can be from the camera before errors render the estimates unreliable. We found that diameter errors do not increase with distance for trees up to 20 m from the camera, but further distances might be possible. Increasing the resolution of the camera will certainly increase the working distance and accuracy of contour localization; however, at the cost of increased processing time. An alternative is to increase the baseline between the cameras; however, this will increase the minimum working distance and trees close to the camera might not appear in both stereo images. Further research is needed to study the relationship between camera configuration, resolution, processing time and accuracy.

6. Patents

US Patent No. US011481972B2: Method of performing dendrometry and forest mapping.

Author Contributions

Conceptualization, L.A.W. and W.C.; methodology, L.A.W. and W.C.; software, L.A.W.; validation, L.A.W. and W.C.; formal analysis, L.A.W.; investigation, L.A.W. and W.C.; resources, L.A.W.; data curation, L.A.W.; writing—original draft preparation, L.A.W.; writing—review and editing, L.A.W. and W.C.; visualization, L.A.W.; supervision, W.C.; project administration, W.C.; funding acquisition, W.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the U.S. Forest Service National Technology and Development Program under contract number 16CS-1113-8100-017.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Lisein, J.; Pierrot-Deseilligny, M.; Bonnet, S.; Lejeune, P. A Photogrammetric Workflow for the Creation of a Forest Canopy Height Model from Small Unmanned Aerial System Imagery. Forests 2013, 4, 922–944. [Google Scholar] [CrossRef]
  2. Suárez, J.C.; Ontiveros, C.; Smith, S.; Snape, S. Use of airborne LiDAR and aerial photography in the estimation of individual tree heights in forestry. Comput. Geosci. 2005, 31, 253–262. [Google Scholar] [CrossRef]
  3. Popescu, S.C.; Wynne, R.H.; Nelson, R.F. Measuring individual tree crown diameter with lidar and assessing its influence on estimating forest volume and biomass. Can. J. Remote Sens. 2003, 29, 564–577. [Google Scholar] [CrossRef]
  4. Erdody, T.L.; Moskal, L.M. Fusion of LiDAR and imagery for estimating forest canopy fuels. Remote Sens. Environ. 2010, 114, 725–737. [Google Scholar] [CrossRef]
  5. Lim, K.; Treitz, P.; Wulder, M.; St-Onge, B.; Flood, M. LiDAR remote sensing of forest structure. Prog. Phys. Geogr. Earth Environ. 2003, 27, 88–106. [Google Scholar] [CrossRef]
  6. Koukoulas, S.; Blackburn, G.A. Mapping individual tree location, height and species in broadleaved deciduous forest using airborne LIDAR and multi-spectral remotely sensed data. Int. J. Remote Sens. 2005, 26, 431–455. [Google Scholar] [CrossRef]
  7. Shelbourne, C.; Namkoong, G. Photogrammetric technique for measuring bole straightness. In Proceedings of the 8th Southern Conference on Forest Tree Improvement, Savannah, GA, USA, 16–17 June 1966; pp. 131–136. [Google Scholar]
  8. Hapca, A.I.; Mothe, F.; Leban, J.M. A digital photographic method for 3D reconstruction of standing tree shape. Ann. For. Sci. 2007, 64, 631–637. [Google Scholar] [CrossRef]
  9. Hapca, A.I.; Mothe, F.; Leban, J.M. Three-dimensional profile classification of standing trees using a stereophotogrammetric method. Scand. J. For. Res. 2008, 23, 46–52. [Google Scholar] [CrossRef]
  10. Côté, J.F.; Fournier, R.A.; Egli, R. An architectural model of trees to estimate forest structural attributes using terrestrial LiDAR. Environ. Model. Softw. 2011, 26, 761–777. [Google Scholar] [CrossRef]
  11. Raumonen, P.; Kaasalainen, M.; Åkerblom, M.; Kaasalainen, S.; Kaartinen, H.; Vastaranta, M.; Holopainen, M.; Disney, M.; Lewis, P. Fast Automatic Precision Tree Models from Terrestrial Laser Scanner Data. Remote Sens. 2013, 5, 491–520. [Google Scholar] [CrossRef]
  12. Maas, H.G.; Bienert, A.; Scheller, S.; Keane, E. Automatic forest inventory parameter determination from terrestrial laser scanner data. Int. J. Remote Sens. 2008, 29, 1579–1593. [Google Scholar] [CrossRef]
  13. Liang, X.; Litkey, P.; Hyyppa, J.; Kaartinen, H.; Vastaranta, M.; Holopainen, M. Automatic Stem Mapping Using Single-Scan Terrestrial Laser Scanning. IEEE Trans. Geosci. Remote Sens. 2012, 50, 661–670. [Google Scholar] [CrossRef]
  14. Olofsson, K.; Holmgren, J. Single Tree Stem Profile Detection Using Terrestrial Laser Scanner Data, Flatness Saliency Features and Curvature Properties. Forests 2016, 7, 207. [Google Scholar] [CrossRef]
  15. Fritz, A.; Kattenborn, T.; Koch, B. UAV-based photogrammetric point clouds – tree stem mapping in open stands in comparison to terrestrial laster scanner point clouds. ISPRS-Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2013, XL-1-W2, 141–146. [Google Scholar] [CrossRef]
  16. Morgenroth, J.; Gomez, C. Assessment of tree structure using a 3D image analysis technique—A proof of concept. Urban For. Urban Green. 2014, 13, 198–203. [Google Scholar] [CrossRef]
  17. Liang, X.; Jaakkola, A.; Wang, Y.; Hyyppä, J.; Honkavaara, E.; Liu, J.; Kaartinen, H. The Use of a Hand-Held Camera for Individual Tree 3D Mapping in Forest Sample Plots. Remote Sens. 2014, 6, 6587–6603. [Google Scholar] [CrossRef]
  18. Miller, J.; Morgenroth, J.; Gomez, C. 3D modelling of individual trees using a handheld camera: Accuracy of height, diameter and volume estimates. Urban For. Urban Green. 2015, 14, 932–940. [Google Scholar] [CrossRef]
  19. Surový, P.; Yoshimoto, A.; Panagiotidis, D. Accuracy of Reconstruction of the Tree Stem Surface Using Terrestrial Close-Range Photogrammetry. Remote Sens. 2016, 8, 123. [Google Scholar] [CrossRef]
  20. Mikita, T.; Janata, P.; Surový, P. Forest Stand Inventory Based on Combined Aerial and Terrestrial Close-Range Photogrammetry. Forests 2016, 7, 165. [Google Scholar] [CrossRef]
  21. Forsman, M.; Börlin, N.; Holmgren, J. Estimation of Tree Stem Attributes Using Terrestrial Photogrammetry with a Camera Rig. Forests 2016, 7, 61. [Google Scholar] [CrossRef]
  22. Fang, R.; Strimbu, B.M. Stem Measurements and Taper Modeling Using Photogrammetric Point Clouds. Remote Sens. 2017, 9, 716. [Google Scholar] [CrossRef]
  23. Kass, M.; Witkin, A.; Terzopoulos, D. Snakes: Active contour models. Int. J. Comput. Vis. 1988, 1, 321–331. [Google Scholar] [CrossRef]
  24. Wells, L.A.; Chung, W. Real-Time Computer Vision for Tree Stem Detection and Tracking. Forests 2023, 14, 267. [Google Scholar] [CrossRef]
  25. He, K.; Gkioxari, G.; Dollár, P.; Girshick, R. Mask R-CNN. arXiv 2018, arXiv:cs.CV/1703.06870. [Google Scholar]
  26. Chung, W.; Lyons, C.K.; Wells, L.A. Innovations in Forest Harvesting Technology. In Chapter 18 In Achieving Sustainable Forestry Volume 1: Boreal and Temperate Forests; Stanturf, J., Ed.; Burleigh Dodds Science Publishing: Cambridge, UK, 2019. [Google Scholar]
  27. Thomas, L.H. Elliptic Problems in Linear Differential Equations over a Network; Technical Report; Columbia University: New York, NY, USA, 1949. [Google Scholar]
  28. Xu, C.; Prince, J.L. Snakes, shapes, and gradient vector flow. IEEE Trans. Image Process. 1998, 7, 359–369. [Google Scholar] [CrossRef]
  29. Hirschmuller, H. Stereo Processing by Semiglobal Matching and Mutual Information. IEEE Trans. Pattern Anal. Mach. Intell. 2008, 30, 328–341. [Google Scholar] [CrossRef]
  30. Hernandez-Juarez, D.; Chacón, A.; Espinosa, A.; Vázquez, D.; Moure, J.C.; López, A.M. Embedded Real-time Stereo Estimation via Semi-Global Matching on the GPU. In Proceedings of the International Conference on Computational Science 2016, ICCS 2016, San Diego, CA, USA, 6–8 June 2016; pp. 143–153. [Google Scholar] [CrossRef]
  31. Wang, Q.; Boyer, K.L. The active geometric shape model: A new robust deformable shape model and its applications. Comput. Vis. Image Underst. 2012, 116, 1178–1194. [Google Scholar] [CrossRef]
  32. Darwin, C. The Power of Movement in Plants; Cambridge University Press: Cambridge, UK, 1880. [Google Scholar] [CrossRef]
  33. Metzger, K. Der Wind als massgebender Faktor für das Wachtsum der Bäume. Mundener Forstl. Hefte 1893, 3, 35–86. [Google Scholar]
  34. Max, T.A.; Burkhart, H.E. Segmented Polynomial Regression Applied to Taper Equations. For. Sci. 1976, 22, 283–289. [Google Scholar]
  35. Kozak, A. A variable-exponent taper equation. Can. J. For. Res. 2011, 18, 1363–1368. [Google Scholar] [CrossRef]
  36. Goodwin, A.N. A cubic tree taper model. Aust. For. 2009, 72, 87–98. [Google Scholar] [CrossRef]
  37. Zhang, Z. A flexible new technique for camera calibration. IEEE Trans. Pattern Anal. Mach. Intell. 2000, 22, 1330–1334. [Google Scholar] [CrossRef]
  38. Fischler, M.A.; Bolles, R.C. Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography. Commun. ACM 1981, 24, 381–395. [Google Scholar] [CrossRef]
  39. Wells, L.A.; Chung, W. Evaluation of Ground Plane Detection for Estimating Breast Height in Stereo Images. For. Sci. 2020, 66, 612–622. [Google Scholar] [CrossRef]
  40. Putra, B.T.W.; Ramadhani, N.J.; Soedibyo, D.W.; Marhaenanto, B.; Indarto, I.; Yualianto, Y. The use of computer vision to estimate tree diameter and circumference in homogeneous and production forests using a non-contact method. For. Sci. Technol. 2021, 17, 32–38. [Google Scholar] [CrossRef]
  41. Ahamed, A.; Foye, J.; Poudel, S.; Trieschman, E.; Fike, J. Measuring Tree Diameter with Photogrammetry Using Mobile Phone Cameras. Forests 2023, 14, 2027. [Google Scholar] [CrossRef]
  42. Eliopoulos, N.J.; Shen, Y.; Nguyen, M.L.; Arora, V.; Zhang, Y.; Shao, G.; Woeste, K.; Lu, Y. Rapid Tree Diameter Computation with Terrestrial Stereoscopic Photogrammetry. J. For. 2020, 118, 355–361. [Google Scholar] [CrossRef]
  43. Shao, T.; Qu, Y.; Du, J. A low-cost integrated sensor for measuring tree diameter at breast height (DBH). Comput. Electron. Agric. 2022, 199, 107140. [Google Scholar] [CrossRef]
  44. NVIDIA Corporation. NVIDIA CUDA C Programming Guide; Version 3.2; NVIDIA: Santa Clara, CA, USA, 2010. [Google Scholar]
  45. Murphy, G. Determining Stand Value and Log Product Yields Using Terrestrial Lidar and Optimal Bucking: A Case Study. J. For. 2008, 106, 317–324. [Google Scholar]
  46. Wang, J.; Liu, J.; LeDoux, C.B. A three-dimensional bucking system for optimal bucking of Central Appalachian hardwoods. Int. J. For. Eng. 2009, 20, 26–35. [Google Scholar] [CrossRef]
  47. Constant, T.; Mothe, F.; Badia, M.A.; Saint-Andre, L. How to relate the standing tree shape to internal wood characteristics: Proposal of an experimental method applied to poplar trees. Ann. For. Sci. 2003, 60, 371–378. [Google Scholar] [CrossRef]
  48. Liu, S.; Bitterlich, W.; Cieszewski, C.J.; Zasada, M.J. Comparing the Use of Three Dendrometers for Measuring Diameters at Breast Height. South. J. Appl. For. 2011, 35, 136–141. [Google Scholar] [CrossRef]
Figure 1. Outline of our proposed method. Given a stereo pair and tree stem detection boxes, the algorithm localizes the edges of the tree in the left and right stereo images and reconstructs their 3D shape. The stems are then mapped to the object coordinate frame [26].
Figure 1. Outline of our proposed method. Given a stereo pair and tree stem detection boxes, the algorithm localizes the edges of the tree in the left and right stereo images and reconstructs their 3D shape. The stems are then mapped to the object coordinate frame [26].
Forests 14 02299 g001
Figure 2. Deformable Stem Model: Five instances of the Deformable Stem Model each with a different number of control points (n). (a) n = 4 ; (b) n = 5 ; (c) n = 6 ; (d) n = 7 ; (e) n = 8 .
Figure 2. Deformable Stem Model: Five instances of the Deformable Stem Model each with a different number of control points (n). (a) n = 4 ; (b) n = 5 ; (c) n = 6 ; (d) n = 7 ; (e) n = 8 .
Forests 14 02299 g002
Figure 3. External energy (images are rescaled for visualization purposes): (a) Gray-scale image of a tree stem, I ( x , y ) . (b) Edge map, I x ( x , y ) . (c) Force map I x x ( x , y ) . (d) Gradient flow force map, G ( x , y ) . (e) Contraction force map, W x ( x , y ) .
Figure 3. External energy (images are rescaled for visualization purposes): (a) Gray-scale image of a tree stem, I ( x , y ) . (b) Edge map, I x ( x , y ) . (c) Force map I x x ( x , y ) . (d) Gradient flow force map, G ( x , y ) . (e) Contraction force map, W x ( x , y ) .
Forests 14 02299 g003
Figure 4. Internal straightness forces: (a) An instance of the Deformable Stem Model exhibiting an unlikely form. (b) Vectors showing the direction and magnitude of the internal straightness forces, S x , scaled up here for visualization purposes. (c) The result of instance (a) after 10 iterations at α = 0.1 with only the internal straightness forces acting on the stem.
Figure 4. Internal straightness forces: (a) An instance of the Deformable Stem Model exhibiting an unlikely form. (b) Vectors showing the direction and magnitude of the internal straightness forces, S x , scaled up here for visualization purposes. (c) The result of instance (a) after 10 iterations at α = 0.1 with only the internal straightness forces acting on the stem.
Forests 14 02299 g004
Figure 5. Internal taper forces: (a) An instance of the Deformable Stem Model exhibiting an unlikely taper function. (b) Vectors showing the direction and magnitude of the internal taper forces, scaled up here for visualization purposes. (c) The result of instance (a) after 10 iterations at β = 0.1 with only the internal taper forces acting on the stem.
Figure 5. Internal taper forces: (a) An instance of the Deformable Stem Model exhibiting an unlikely taper function. (b) Vectors showing the direction and magnitude of the internal taper forces, scaled up here for visualization purposes. (c) The result of instance (a) after 10 iterations at β = 0.1 with only the internal taper forces acting on the stem.
Forests 14 02299 g005
Figure 6. Diagram showing the flow of operations in the Deformable Stem Model algorithm.
Figure 6. Diagram showing the flow of operations in the Deformable Stem Model algorithm.
Forests 14 02299 g006
Figure 7. Six examples of optimization of the Deformable Stem Model on real images. (a) Homogeneously illuminated stem with clear edges. (b) Discontinuous edges due to branching. (c) Harsh shadows crossing the stem. (d) Fully illuminated stem with low contrast to background. (e) Partially occluded stem in the background. (f) Curved stem at maximum working distance of camera. (g) Partially illuminated stem. (h) Low contrast stem with branching.
Figure 7. Six examples of optimization of the Deformable Stem Model on real images. (a) Homogeneously illuminated stem with clear edges. (b) Discontinuous edges due to branching. (c) Harsh shadows crossing the stem. (d) Fully illuminated stem with low contrast to background. (e) Partially occluded stem in the background. (f) Curved stem at maximum working distance of camera. (g) Partially illuminated stem. (h) Low contrast stem with branching.
Forests 14 02299 g007
Figure 8. Diameter estimation results. (a): Diameter vs. projected width, 2 τ ( y ) , at observed breast height, y, scaled by observed distance, z. The gray dashed line is the best-fit linear model. (b): Observed distance, z, vs. predicted distance at observed breast height, z ^ ( y ) , and predicted distance at predicted breast height, z ^ ( y ^ ) . (c): Observed diameter, d, vs. predicted diameter at predicted distance and predicted breast height, d ^ ( z ^ , y ^ ) . (d): Observed distance, z vs. diameter residuals, d ^ ( z ^ , y ^ ) d .
Figure 8. Diameter estimation results. (a): Diameter vs. projected width, 2 τ ( y ) , at observed breast height, y, scaled by observed distance, z. The gray dashed line is the best-fit linear model. (b): Observed distance, z, vs. predicted distance at observed breast height, z ^ ( y ) , and predicted distance at predicted breast height, z ^ ( y ^ ) . (c): Observed diameter, d, vs. predicted diameter at predicted distance and predicted breast height, d ^ ( z ^ , y ^ ) . (d): Observed distance, z vs. diameter residuals, d ^ ( z ^ , y ^ ) d .
Forests 14 02299 g008
Table 1. Univariate statistics of diameter observations. Units are in centimeters.
Table 1. Univariate statistics of diameter observations. Units are in centimeters.
ObservationMeanSTDMinMax
Distance (z)1017.45458.85350.521813.56
Breast ht. (y)137.160.00137.16137.16
Diameter (d)35.094.8826.0042.30
Table 2. Deformable Stem Model hyperparameters used in validation experiment.
Table 2. Deformable Stem Model hyperparameters used in validation experiment.
Hyper-ParameterDescriptionValue
nnumber of control points5
α internal straightness force step size0.1
β internal taper force step size0.1
κ contraction force scaling parameter0.1
μ stereo smoothness constraint parameter0.3
Table 3. RMSE, bias, and MAE for the distance, breast height, and diameter estimators. The bold values show the results for the fully automatic diameter estimator.
Table 3. RMSE, bias, and MAE for the distance, breast height, and diameter estimators. The bold values show the results for the fully automatic diameter estimator.
EstimatorEquationMeanSTDRMSEBiasMAE
(cm)(cm)(cm)(cm)(cm)
Distance z ^ ( y ) Equation (47)1023.34468.3833.005.8923.26
z ^ ( y ^ ) Equation (47)1023.82468.4733.276.3823.10
Breast ht. y ^ Equation (49)135.035.856.10 2.13 4.87
Diameter d ^ ( z , y ) Equation (46)35.115.231.610.021.39
d ^ ( z ^ , y ) Equation (48)35.225.511.760.131.49
d ^ ( z ^ , y ^ ) Equation (50)35.395.421.940.301.55
Table 4. Run time results for various components of the algorithm.
Table 4. Run time results for various components of the algorithm.
Size (kp)0–1010–2020–3030–4040–50
Upload (ms)0.040.040.070.080.10
Allocate (ms)0.160.240.400.580.77
Gradient flow (ms)0.881.121.912.703.57
Contraction force (ms)0.240.320.440.580.70
Optimization-1002.423.204.184.775.18
iterations (ms)
Total (ms)3.744.936.998.7110.30
Frames per second267.34202.85142.99114.8097.05
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Wells, L.A.; Chung, W. A Deformable Shape Model for Automatic and Real-Time Dendrometry. Forests 2023, 14, 2299. https://doi.org/10.3390/f14122299

AMA Style

Wells LA, Chung W. A Deformable Shape Model for Automatic and Real-Time Dendrometry. Forests. 2023; 14(12):2299. https://doi.org/10.3390/f14122299

Chicago/Turabian Style

Wells, Lucas A., and Woodam Chung. 2023. "A Deformable Shape Model for Automatic and Real-Time Dendrometry" Forests 14, no. 12: 2299. https://doi.org/10.3390/f14122299

APA Style

Wells, L. A., & Chung, W. (2023). A Deformable Shape Model for Automatic and Real-Time Dendrometry. Forests, 14(12), 2299. https://doi.org/10.3390/f14122299

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop