Wood Quality Assessment of Standing Stems When Measurements Do Not Line Up: A Knot Geometry Model for Spatially Misaligned Tree Data

1. Introduction

In North America, the assessment of standing wood quality has emerged as a critical early intervention strategy in modern forest management, marking a fundamental shift from traditional volume-based metrics to value-based recovery [1,2]. This paradigm shift is essential for optimizing the wood supply chain, ensuring the structural soundness of timber products, and maximizing market value while simultaneously reducing waste and supporting sustainable silvicultural practices [3,4,5]. By implementing assessment protocols early in the production cycle, forest managers can facilitate the strategic segregation of timber by quality prior to harvest. For example, adding only simple visual parameters like number of branches in the lower part of the bole has significantly improved prediction of lumber yield [6,7]. This proactive approach enables the precise allocation of raw materials to their most economically viable end-uses, such as distinguishing high-stiffness structural lumber from pulpwood, thereby maximizing the return on investment for each individual stem [3,2].

Beyond immediate economic gains, early assessment is a prerequisite for meeting the rigorous mechanical standards for strength and stiffness required by the modern construction industry [1,4]. Instead of viewing timber as a generic product, managers must align forest growth with industry demands by predicting wood quality long before harvest. Predicting the quality of potential products within specific growing environments allows managers to tailor silvicultural prescriptions to optimize both the quantity and internal characteristics of the standing resource [8,9]. Because the standing tree represents the origin of the forest products value chain, it is the point at which wood quality is most susceptible to management influence [10]. Consequently, targeting wood quality during silvicultural planning creates value that extends beyond mere fiber extraction, allowing for more robust decision-making under varying environmental and management scenarios [3,8]. However, to realize these benefits across broad forest landscapes, it is necessary to move beyond plot-level sampling toward integrated, large-scale assessment methodologies

In recent years, LiDAR (Light Detection and Range) has evolved from a simple inventory tool and volume assessment to a sophisticated technology that either substitutes or complements traditional Non-Destructive Evaluation (NDE) by predicting internal traits via external structure and high-resolution geometry [2]. Although current resolutions are insufficient to detect individual branch scars, high resolution LiDAR have been successfully employed to delineate crown features, canopy density, bole characteristics such as crown base height (CBH), diameter along the bole, stem straightness (sweep) etc. These metrics serve as effective proxies for wood quality in various stand types [11,12,13,14]. Furthermore, airborne LiDAR (ALS) is seen to be effective for predicting site-level wood density and general stem dimensions across large territories [14]. However, ALS often suffers from significant occlusion near-nadir scanning angles and canopy interference, which prevents detailed characterization of the lower stem where the most valuable sawlogs are located [15]. Operating at lower altitudes, slower speed and multiple scanning perspectives, drone or UAV based LiDAR ( ULS), with its significantly high density can penetrate the canopy more effectively to capture lower-stem morphology [13,16].

The size and distribution of knots are widely recognized as critical determinants of wood quality, directly influencing lumber grading, mechanical performance, and aesthetic value [2,17,18]. Although restricted to plot-level assessments, the proximity and perspective of Terrestrial (TLS) and Mobile Laser Scanning (MLS) enable the capture of millimeter-level point clouds. These high-density datasets allow for the accurate 3D reconstruction of woody structures, including the measurement of taper, sweep, and ovality [19,20], as well as the detection of branch scars and whorl patterns [21,22].

Recent studies have leveraged these external morphological features to infer internal characteristics; for instance, TLS data has been correlated with CT scans to achieve 62.5% accuracy in the automated quality grading of European beech sawlogs [23]. Furthermore, when integrated with resistance drilling, these scanners can predict wood basic density with moderate success of an R² of 0.51 [14]. While they provide high-resolution knot-level detail required for internal trait prediction, their limited spatial reach necessitates a bridge toward more scalable platforms for operational efficiency. One strategy for scaling wood quality predictions could be by leveraging ULS-derived crown and bole metrics as inputs for empirical models. These models, calibrated using meticulous ground-based data, allow for the translation of external structural geometry into broad-scale wood quality indicators.

Despite the high-resolution capabilities, spatial misalignment between external predictors and internal responses remains a fundamental bottleneck when integrating multi-source datasets [24]. In wood quality modeling, the covariates derived from LiDAR (such as external bole surface irregularities or branch metrics along the bole) and the response variables (such as internal knot geometry from destructive sampling or Computed Tomography - CT scans) are rarely evaluated at the same spatial coordinates. Despite efforts, this “change of support" problem introduces significant bias, as traditional regression techniques often assume a one-to-one spatial correspondence between variables [25,26,27]. Standard alignment methods frequently fail to account for the irregular 3D topology of the tree or the sensor-specific occlusions that shift the perceived location of a feature [28,29]. Consequently, there is a critical need for a generalized modeling framework capable of reconciling these disjointed data points without sacrificing the localized precision required for knot-level analysis or even log-segment level analysis.

The primary objective of this study is to develop a generalized modeling framework for cases where the covariate (e.g., external tree characteristics) and the response variable (e.g., knot geometry) are evaluated at disjointed spatial points. By leveraging the Law of Total Expectation in conjunction with Kernel Density Estimation, we propose a novel statistical method, called regression of misaligned covariates (RMC), that reconciles these misaligned datasets, enabling the prediction of internal wood quality attributes from external structural tree geometry.

The development of such a framework requires validation across contrasting physiological profiles; thus, red pine (Pinus resinosa) and eastern white pine (Pinus strobus) are utilized here as distinct prototypes of conifer physiology. Red pine, characterized by exceptional genetic and structural uniformity, serves as a stable baseline for quality-related assessments [30,31]. In contrast, eastern white pine presents significant geometric complexity due to its discrete branching architecture and susceptibility to terminal leader damage [32,33,34]. These biological stressors often result in morphological irregularities, including crooks, sweeps, and multi-leader crowns as well as positioning the species as a rigorous testbed for LiDAR-based algorithms designed to map knot distribution and volumetric recovery [35,36].

Through this dual-species validation, we illustrate the estimation of primary wood quality metrics in standing trees and demonstrate how these data can directly inform management decisions and industrial processing. By bridging the gap between raw spatial data and volumetric recovery, this research ultimately provides a robust, scalable pathway for assessing wood quality at both the individual tree and sub-tree levels across extensive forest blocks.

2. Materials and Methods

2.1. Regression of Misaligned Covariates (RMC)

2.1.1. Model Presentation

The goal is to predict the knot geometry, such as volume, diameter and size, using tree external characteristics such as log diameter, branch length and branch diameter measured along the height (or length) of the tree bole. However, in contrast with the usual regression framework, here we deal with misaligned covariates i.e., when the data for the covariates and the response variables are collected at different times and at varying heights along the bole. Specifically:

• the knot geometry measurements (internal characteristics) are available at certain height points.
• the external tree characteristics are recorded at other, possibly different, height points.

Given Y as a response variable (indexed by 0) as tree external characteristics, X as a covariate (indexed by 1) as external characteristics, and H as the relative height (indexed by 2), the sample measurements could be denoted as:

$${\left\{\left(Y_i,H_i\right)\right\}}_{i=1}^{i=n_{02}} \mathrm{and} {\left\{\left(X_i,H_i\right)\right\}}_{i=1}^{i=n_{12}}$$

The objective, hence, is to $\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{k}Y$ to $X$, taking note that the $Y_i$ 's to the $X_i$ 's are misaligned i.e. are not available at the same height points. Based on the available data, one can then estimate

$$m_{02}\left(h\right) :=\mathbb{E}(Y\mid H=h)$$

(1)

$$m_{12}\left(h\right) :=\mathbb{E}(X\mid H=h)$$

(2)

The goal is to estimate its regression given as:

$$m_{012}(x,h):=\mathbb{E}(Y\mid X=x,H=h)$$

(3)

The proposed approach here is based on the well-known result that the expectation of conditional expectation is an unconditional expectation, or the law of total expectation, and given by the following equations:

$$\{\begin{array}{}\mathrm{(4)}& {\mathbb{E}}_X\left[\mathbb{E}\right(Y\mid X,h)\mid H=h]=\mathbb{E}[Y\mid H=h]\mathrm{(5)}& {\mathbb{E}}_{X, H}\left[\mathbb{E}\right(Y\mid X,H\left)\right]=\mathbb{E}\left[Y\right] \end{array}$$

For its estimation, a suitable finite-dimensional model for $m_{012}\left(x,h\right)$ could be chosen, for example a form like:

$$m_{012}(x,h)=a_0+a_{1}x+a_{2}h+a_3{x}^2+a_4{h}^2+\dots$$

(6)

and solve for the unknowns $a_0,a_1,\dots$Equations (4) and (5) can be written more precisely as,

$$\{\begin{array}{}\mathrm{(7)}& \int \left(\int {\displaystyle \frac{y{f}_{012}(y,x,h)\mathrm{d}y}{f_{12}(x,h)}}\right){\displaystyle \frac{f_{12}(x,h)\mathrm{d}x}{f_2\left(h\right)}}=\int {\displaystyle \frac{y{f}_{02}(y,h)\mathrm{d}y}{f_2\left(h\right)}}\mathrm{(8)}& \iint \left(\int {\displaystyle \frac{y{f}_{012}(y,x,h)\mathrm{d}y}{f_{12}(x,h)}}\right)f_{12}(x,h)\mathrm{d}x \mathrm{d}h=\iiint y{f}_{012}(y,x,h)\mathrm{d}y \mathrm{d}x \mathrm{d}h\end{array}$$

where it is assumed that

$(Y,X,H)$ has joint density $f_{012}(y,x,h)$$(Y,H)$ has joint density $f_{02}(y,h)$$(X,H)$ has joint density $f_{12}(x,h)$$H$ has marginal density $f_2\left(h\right)$,

and

$$\int f_{012}(y,x,h)\mathrm{d}x=f_{02}(y,h)$$

(9)

$$\mathrm{and} \int f_{012}(y,x,h)\mathrm{d}x \mathrm{d}h=f_0\left(y\right)$$

(10)

$f_0\left(y\right)$ being the marginal density of $Y$.

2.1.2. Regression Estimation Approach

Cancelling $f_2\left(h\right)$ from Equation (7), one can write

$$\int \left(\underset{m_{012}(x,h)}{\underbrace{\int {\displaystyle \frac{y{f}_{012}(y,x,h)\mathrm{d}y}{f_{12}(x,h)}}}}\right){\displaystyle \frac{f_{12}(x,h)\mathrm{d}x}{f_2\left(h\right)}}=\int {\displaystyle \frac{y{f}_{02}(y,h)\mathrm{d}y}{f_2\left(h\right)}} \mathrm{for} \mathrm{all} h$$

(11)

$$⟺\int m_{012}(x,h)f_{12}(x,h)\mathrm{d}x=\int y{f}_{02}(y,h)\mathrm{d}y \mathrm{for} \mathrm{all} h$$

(12)

so that Equations (7) and (8) become

$$\{\begin{array}{}\mathrm{(13)}& \int m_{012}\left(x,h\right)f_{12}\left(x,h\right)dx=\int y{f}_{02}\left(y,h\right)dy\mathrm{(14)}& \mathbb{E}\left[m_{012}(X,H)\right]=\mathbb{E}\left(Y\right)\end{array}$$

Plugging in the corresponding kernel estimators on both sides we have

$$\{\begin{array}{}\mathrm{(15)}& \int m_{012}\left(x,h\right){\hat{f}}_{12}\left(x,h\right)dx=\int y{\hat{f}}_{02}\left(y,h\right)dy\mathrm{(16)}& {\displaystyle \frac{1}{n_{12}}}\sum _{i=1}^{n_{12}} m_{012}\left(X_i,H_i\right)={\displaystyle \frac{1}{n_{02}}}\sum _{i=1}^{n_{02}} Y_i\end{array}$$

where

$$\int {\hat{f}}_{12}(x,h)\mathrm{d}x={\displaystyle \frac{1}{n_{12}{b}_n^2}}\sum _{i=1}^{n_{12}} K\left({\displaystyle \frac{x-X_i}{b_n}}\right)K\left({\displaystyle \frac{h-H_i}{b_n}}\right) \mathrm{for} j=1,\dots ,n_{12}$$

(17)

$$\int {\hat{f}}_{02}(y,h)\mathrm{d}x={\displaystyle \frac{1}{n_{02}{b}_n^2}}\sum _{i=1}^{n_{02}} K\left({\displaystyle \frac{y-Y_i}{b_n}}\right)K\left({\displaystyle \frac{h-H_i}{b_n}}\right) \mathrm{for} j=1,\dots ,n_{02}$$

(18)

and $K($.$)is a kernel.$ See, for example, [37] for an introduction to kernel smoothing.

The above equations (15) and (16) can be solved by choosing a suitable finite-dimensional model for

$$m_{012}\left(x,h\right)$$

as mentioned above, for example, assuming

$$m(x,h)=a_0+a_{1}x+a_{2}h+a_3{x}^2+a_4{h}^2+\dots$$

(19)

Coefficients $a_0,a_1,\dots$, can be solved by creating as many equations as the number of coefficients by considering different values of $h:h_0,h_1,h_2,\dots$To illustrate the proposed estimation method, the following subsection examines a tri-variate normal distribution of the variables of interest.

2.1.3. (Y,X,H) Follow a 3- Variate Normal

Suppose

$$\left(\begin{array}{l}Y\\ X\\ H\end{array}\right)\sim {\mathcal{N}}_3\left(\left(\begin{array}{l}{\mu }_0\\ {\mu }_1\\ {\mu }_2\end{array}\right),\mathsf{\Sigma }\right)$$

(20)

From the properties of multivariate Normal distributions, the expectations $m_{12}\left(h\right)=\mathbb{E}(X\mid H=h),m_{02}\left(h\right)=\mathbb{E}(Y\mid H=h)$ and $m(x,h)=m_{012}(x,h)=\mathbb{E}(Y\mid X=x,H=h)$ are linear.

Thus we may write

$$m_{12}\left(h\right)=\mathbb{E}(X\mid H=h)=b_0+b_{2}h$$

(21)

$$\mathrm{and} m_{02}\left(h\right)=\mathbb{E}(Y\mid H=h)=c_0+c_{2}h.$$

(22)

And

$$m(x,h)=m_{012}(x,h)=\mathbb{E}(Y\mid X=x,H=h)=a_0+a_{1}x+a_{2}h$$

(23)

Equations (4) and (5) then become

$$\{\begin{array}{}\mathrm{(24)}& a_0+a_1\left(b_0+b_{2}h\right)+a_{2}h=c_0+c_{2}h, \mathrm{for} \mathrm{all} h\mathrm{(25)}& a_0+a_1{\mu }_1+a_2{\mu }_2={\mu }_0\end{array}$$

This yields a system of three equations with three unknowns ( $a_0,a_1,a_2$ ):

$$\{\begin{array}{}\mathrm{(26)}& a_0+b_0{a}_1=c_0\mathrm{(27)}& b_2{a}_1+a_2=c_2\mathrm{(28)}& a_0+a_1{\mu }_1+a_2{\mu }_2={\mu }_0\end{array}$$

That can be represented in a matrix form as:

$$\left(\begin{array}{ccc}1& b_0& 0\\ 0& b_2& 1\\ 1& {\mu }_1& {\mu }_2\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right)=\left(\begin{array}{c}{c}_0\\ c_2\\ {\mu }_0\end{array}\right)$$

(29)

$$\left(\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right)={\left(\begin{array}{ccc}1& b_0& 0\\ 0& b_2& 1\\ 1& {\mu }_1& {\mu }_2\end{array}\right)}^{-1}\left(\begin{array}{c}{c}_0\\ c_2\\ {\mu }_0\end{array}\right)$$

(30)

Given a sample,

$${\hat{\mu }}_0={\displaystyle \frac{1}{n_{02}}}\sum _{i=1}^{n_{02}} Y_i$$

(31)

$${\hat{\mu }}_1={\displaystyle \frac{1}{n_{12}}}\sum _{i=1}^{n_{12}} X_i$$

(32)

$${\hat{\mu }}_2={\displaystyle \frac{1}{n_2}}\sum _{i=1}^{n_2} H_i={\displaystyle \frac{1}{n_{02}+n_{12}}}\left(\sum _{i=1}^{n_{02}} H_i+\sum _{i=1}^{n_{12}} H_i\right)$$

(33)

where $b_0,b_2,c_0$ and $c_2$ can be estimated from Equations (21) and (22) using Ordinary Least Squares (OLS) technique. This should yield

$$\left(\begin{array}{l}{a}_0\\ a_1\\ a_2\end{array}\right)={\left(\begin{array}{ccc}1& {\hat{b}}_0& 0\\ 0& {\hat{b}}_2& 1\\ 1& {\hat{\mu }}_1& {\hat{\mu }}_2\end{array}\right)}^{-1}\left(\begin{array}{c}{\hat{c}}_0\\ {\hat{c}}_2\\ {\hat{\mu }}_0\end{array}\right)$$

(34)

2.1.4. Generalisation of the Model

To demonstrate the flexibility of the proposed framework, three specific functional forms for m (x,h) are evaluated and their estimations presented in this section.

Case 1: m (x,h) is linear

Given m (x,h) is linear,

$$m_{012}(x,h)=\mathbb{E}(Y\mid X=x,H=h)=a_0+a_{1}x+a_{2}h$$

(35)

And equations (15) and (16) become:

$$\{\begin{array}{}\mathrm{(36)}& \int \left(a_0+a_{1}x+a_{2}h\right){\hat{f}}_{12}\left(x,h\right)\mathrm{d}x=\int y{\hat{f}}_{02}\left(y,h\right)\mathrm{d}y\mathrm{(37)}& {\displaystyle \frac{1}{n_{12}}}\sum _{i=1}^{n_{12}} \left(a_0+a_1{X}_i+a_2{H}_i\right)={\displaystyle \frac{1}{n_{02}}}\sum _{i=1}^{n_{02}} Y_i\end{array}$$

$$\{\begin{array}{}\mathrm{(38)}& \int \left(a_0+a_{1}x+a_{2}h\right){\hat{f}}_{12}\left(x,h\right)\mathrm{d}x=\int y{\hat{f}}_{02}\left(y,h\right)\mathrm{d}y\mathrm{(39)}& a_0+a_1{\bar{X}}_{n_{12}}+a_2{\bar{H}}_{n_{12}}={\bar{Y}}_{n_{02}}\end{array}$$

where

$${\bar{X}}_{n_{12}}=\iint x{\hat{f}}_{12}(x,h)\mathrm{d}x \mathrm{d}h={\displaystyle \frac{1}{n_{12}}}\sum _{i=1}^{n_{12}} X_i$$

(40)

$${\bar{H}}_{n_{12}}=\iint h{\hat{f}}_{12}(x,h)\mathrm{d}x \mathrm{d}h={\displaystyle \frac{1}{n_{12}}}\sum _{i=1}^{n_{12}} H_i$$

(41)

$${\bar{Y}}_{n_0}=\int y{\hat{f}}_2\left(y\right)\mathrm{d}y={\displaystyle \frac{1}{n_0}}\sum _{i=1}^{n_0} Y_i$$

(42)

Taking $h=h_1$ and $h=h_2$in the first equation (Equation 36), we have

$$\{\begin{array}{}\mathrm{(43)}& \int \left(a_0+a_{1}x+a_2{h}_1\right){\hat{f}}_{12}\left(x,h_1\right)\mathrm{d}x=\int y{\hat{f}}_{02}\left(y,h_1\right)\mathrm{d}y\mathrm{(44)}& \int \left(a_0+a_{1}x+a_2{h}_2\right){\hat{f}}_{12}\left(x,h_2\right)\mathrm{d}x=\int y{\hat{f}}_{02}\left(y,h_2\right)\mathrm{d}y\mathrm{(45)}& a_0+a_1{\bar{X}}_{n_{12}}+a_2{\bar{H}}_{n_{12}}={\bar{Y}}_{n_{02}}\end{array}$$

$$\{\begin{array}{}\mathrm{(46)}& a_0\int {\hat{f}}_{12}\left(x,h_1\right)\mathrm{d}x+a_1\int x{\hat{f}}_{12}\left(x,h_1\right)\mathrm{d}x+a_2{h}_1\int {\hat{f}}_{12}\left(x,h_1\right)\mathrm{d}x=\int y{\hat{f}}_{02}\left(y,h_1\right)\mathrm{d}y\mathrm{(47)}& a_0\int {\hat{f}}_{12}\left(x,h_2\right)\mathrm{d}x+a_1\int x{\hat{f}}_{12}\left(x,h_2\right)\mathrm{d}x+a_2{h}_1\int {\hat{f}}_{12}\left(x,h_2\right)\mathrm{d}x=\int y{\hat{f}}_{02}\left(y,h_2\right)\mathrm{d}y\mathrm{(48)}& a_0+a_1{\bar{X}}_{n_{12}}+a_2{\bar{H}}_{n_{12}}={\bar{Y}}_{n_{02}}\end{array}$$

$$\{\begin{array}{}\mathrm{(49)}& a_0{\hat{f}}_2\left(h_1\right)+a_1\int x{\hat{f}}_{12}\left(x,h_1\right)\mathrm{d}x+a_2{h}_1{\hat{f}}_2\left(h_1\right)=\int y{\hat{f}}_{02}\left(y,h_1\right)\mathrm{d}y\mathrm{(50)}& a_0{\hat{f}}_2\left(h_2\right)+a_1\int x{\hat{f}}_{12}\left(x,h_2\right)\mathrm{d}x+a_2{h}_2{\hat{f}}_2\left(h_2\right)=\int y{\hat{f}}_{02}\left(y,h_2\right)\mathrm{d}y\mathrm{(51)}& a_0+a_1{\bar{X}}_{n_{12}}+a_2{\bar{H}}_{n_{12}}={\bar{Y}}_{n_{02}}\end{array}$$

$$\left(\begin{array}{ccc}{\hat{f}}_2\left(h_1\right)& \int x{\hat{f}}_{12}\left(x,h_1\right)\mathrm{d}x& h_1{\hat{f}}_2\left(h_1\right)\\ {\hat{f}}_2\left(h_2\right)& \int x{\hat{f}}_{12}\left(x,h_2\right)\mathrm{d}x& h_2{\hat{f}}_2\left(h_2\right)\\ 1& & {\bar{X}}_{n_{12}}\\ & & {\bar{H}}_{n_{12}}\end{array}\right)\left(\begin{array}{l}{a}_0\\ a_1\\ a_2\end{array}\right)=\left(\begin{array}{c}\int y{\hat{f}}_{02}\left(y,h_1\right)\mathrm{d}y\\ \int y{\hat{f}}_{02}\left(y,h_2\right)\mathrm{d}y\\ {\bar{Y}}_{n_{02}}\end{array}\right)$$

(52)

$$\left(\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right)={\left(\begin{array}{ccc}{\hat{f}}_2\left(h_1\right)& \int x{\hat{f}}_{12}\left(x,h_1\right)\mathrm{d}x& h_1{\hat{f}}_2\left(h_1\right)\\ {\hat{f}}_2\left(h_2\right)& \int x{\hat{f}}_{12}\left(x,h_2\right)\mathrm{d}x& h_2{\hat{f}}_2\left(h_2\right)\\ 1& {\bar{X}}_{n_{12}}& {\bar{H}}_{n_{12}}\end{array}\right)}^{-1}\left(\begin{array}{c}\int y{\hat{f}}_{02}\left(y,h_1\right)\mathrm{d}y\\ \int y{\hat{f}}_{02}\left(y,h_2\right)\mathrm{d}y\\ {\bar{Y}}_{n_{02}}\end{array}\right)$$

(53)

where

$${\hat{f}}_2\left(h_j\right)={\displaystyle \frac{1}{n_2{b}_n}}\sum _{i=1}^{n_2} K\left({\displaystyle \frac{h_j-H_i}{b_n}}\right) \mathrm{for} j=1,\dots ,n_2$$

(54)

$$\int x{\hat{f}}_{12}\left(x,h_j\right)\mathrm{d}x={\displaystyle \frac{1}{n_{12}{b}_n}}\sum _{i=1}^{n_{12}} X_{i}K\left({\displaystyle \frac{h_j-H_i}{b_n}}\right) \mathrm{for} j=1,\dots ,n_{12}$$

(55)

$$\int y{\hat{f}}_{02}\left(y,h_j\right)\mathrm{d}y={\displaystyle \frac{1}{n_{02}{b}_n}}\sum _{i=1}^{n_{02}} Y_{i}K\left({\displaystyle \frac{h_j-H_i}{b_n}}\right) \mathrm{for} j=1,\dots ,n_{02}$$

(56)

Case 2: m (x,h) is quadratic

When m (x,h) is quadratic,

$$m_{012}(x,h)=\mathbb{E}(Y\mid X=x,H=h)=a_0+a_{1}x+a_{2}h+a_3{x}^2+a_4{h}^2+a_{5}xh$$

(57)

And equations (15) and (16) become:

$$\{\begin{array}{}\mathrm{(58)}& \int \left(a_0+a_{1}x+a_{2}h+a_3{x}^2+a_4{h}^2+a_{5}xh\right){\hat{f}}_{12}\left(x,h\right)\mathrm{d}x=\int y{\hat{f}}_{02}\left(y,h\right)\mathrm{d}y\mathrm{(59)}& {\displaystyle \frac{1}{n_{12}}}\sum _{i=1}^{n_{12}} \left(a_0+a_1{X}_i+a_2{H}_i+a_3{X}_i^2+a_4{H}_i^2+a_5{X}_i{H}_i\right)={\displaystyle \frac{1}{n_{02}}}\sum _{i=1}^{n_{02}} Y_i\end{array}$$

(58)

$$\{\begin{array}{}\mathrm{(60)}& a_0{\hat{f}}_2\left(h\right)+a_1\int x{\hat{f}}_{12}(x,h)\mathrm{d}x+a_{2}h{\hat{f}}_2\left(h\right)+a_3\int x^2{\hat{f}}_{12}(x,h)\mathrm{d}x+a_4{h}^2{\hat{f}}_2\left(h\right)+a_{5}h\int x{\hat{f}}_{12}(x,h)\mathrm{d}x=\int y{\hat{f}}_{02}(y,h)\mathrm{d}y\mathrm{(61)}& a_0+a_1{\bar{X}}_{n_{12}}+a_2{\bar{H}}_{n_{12}}+{\displaystyle \frac{a_3}{n_{12}}}\sum _{i=1}^{n_{12}} X_i^2+{\displaystyle \frac{a_4}{n_{12}}}\sum _{i=1}^{n_{12}} H_i^2+{\displaystyle \frac{a_5}{n_{12}}}\sum _{i=1}^{n_{12}} X_i{H}_i={\bar{Y}}_{n_{02}}\end{array}$$

Taking $h=h_1$ and $h=h_2$ for the first equation, we have

$$\{\begin{array}{}\mathrm{(62)}& a_0{\hat{f}}_2\left(h_1\right)+a_1\int x{\hat{f}}_{12}\left(x,h_1\right)\mathrm{d}x+a_2{h}_1{\hat{f}}_2\left(h_1\right)+a_3\int x^2{\hat{f}}_{12}\left(x,h_1\right)\mathrm{d}x+a_4{h}_1^2{\hat{f}}_2\left(h_1\right)+a_5{h}_1\int x{\hat{f}}_{12}\left(x,h_1\right)\mathrm{d}x=\int y{\hat{f}}_{02}\left(y,h_1\right)\mathrm{d}y\mathrm{(63)}& a_0{\hat{f}}_2\left(h_2\right)+a_1\int x{\hat{f}}_{12}\left(x,h_2\right)\mathrm{d}x+a_2{h}_2{\hat{f}}_2\left(h_2\right)+a_3\int x^2{\hat{f}}_{12}\left(x,h_2\right)\mathrm{d}x+a_4{h}_2^2{\hat{f}}_2\left(h_2\right)+a_5{h}_2\int x{\hat{f}}_{12}\left(x,h_2\right)\mathrm{d}x=\int y{\hat{f}}_{02}\left(y,h_2\right)\mathrm{d}y\mathrm{(64)}& a_0{\hat{f}}_2\left(h_3\right)+a_1\int x{\hat{f}}_{12}\left(x,h_3\right)\mathrm{d}x+a_2{h}_3{\hat{f}}_2\left(h_3\right)+a_3\int x^2{\hat{f}}_{12}\left(x,h_3\right)\mathrm{d}x+a_4{h}_3^2{\hat{f}}_2\left(h_3\right)+a_5{h}_3\int x{\hat{f}}_{12}\left(x,h_3\right)\mathrm{d}x=\int y{\hat{f}}_{02}\left(y,h_3\right)\mathrm{d}y\mathrm{(65)}& a_0{\hat{f}}_2\left(h_4\right)+a_1\int x{\hat{f}}_{12}\left(x,h_4\right)\mathrm{d}x+a_2{h}_4{\hat{f}}_2\left(h_4\right)+a_3\int x^2{\hat{f}}_{12}\left(x,h_4\right)\mathrm{d}x+a_4{h}_4^2{\hat{f}}_2\left(h_4\right)+a_5{h}_4\int x{\hat{f}}_{12}\left(x,h_4\right)\mathrm{d}x=\int y{\hat{f}}_{02}\left(y,h_4\right)\mathrm{d}y\mathrm{(66)}& a_0{\hat{f}}_2\left(h_5\right)+a_1\int x{\hat{f}}_{12}\left(x,h_5\right)\mathrm{d}x+a_2{h}_5{\hat{f}}_2\left(h_5\right)+a_3\int x^2{\hat{f}}_{12}\left(x,h_5\right)\mathrm{d}x+a_4{h}_5^2{\hat{f}}_2\left(h_5\right)+a_5{h}_5\int x{\hat{f}}_{12}\left(x,h_5\right)\mathrm{d}x=\int y{\hat{f}}_{02}\left(y,h_5\right)\mathrm{d}y\mathrm{(67)}& a_0+a_1{\bar{X}}_{n_{12}}+a_2{\bar{H}}_{n_{12}}+{\displaystyle \frac{a_3}{n_{12}}}\sum _{i=1}^{n_{12}} X_i^2+{\displaystyle \frac{a_4}{n_{12}}}\sum _{i=1}^{n_{12}} H_i^2+{\displaystyle \frac{a_5}{n_{12}}}\sum _{i=1}^{n_{12}} X_i{H}_i={\bar{Y}}_{n_{02}}\end{array}$$

In matrix form

$$\left(\begin{array}{cccccc}{\hat{f}}_2\left(h_1\right)& \int x{\hat{f}}_{12}\left(x,h_1\right)\mathrm{d}x& h_1{\hat{f}}_2\left(h_1\right)& \int x^2{\hat{f}}_{12}\left(x,h_1\right)\mathrm{d}x& h_1^2{\hat{f}}_2\left(h_1\right)& h_1\int x{\hat{f}}_{12}\left(x,h_1\right)\mathrm{d}x\\ {\hat{f}}_2\left(h_2\right)& \int x{\hat{f}}_{12}\left(x,h_2\right)\mathrm{d}x& h_2{\hat{f}}_2\left(h_2\right)& \int x^2{\hat{f}}_{12}\left(x,h_2\right)\mathrm{d}x& h_2^2{\hat{f}}_2\left(h_2\right)& h_2\int x{\hat{f}}_{12}\left(x,h_2\right)\mathrm{d}x\\ {\hat{f}}_2\left(h_3\right)& \int x{\hat{f}}_{12}\left(x,h_3\right)\mathrm{d}x& h_3{\hat{f}}_2\left(h_3\right)& \int x^2{\hat{f}}_{12}\left(x,h_3\right)\mathrm{d}x& h_3^2{\hat{f}}_2\left(h_3\right)& h_3\int x{\hat{f}}_{12}\left(x,h_3\right)\mathrm{d}x\\ {\hat{f}}_2\left(h_4\right)& \int x{\hat{f}}_{12}\left(x,h_4\right)\mathrm{d}x& h_4{\hat{f}}_2\left(h_4\right)& \int x^2{\hat{f}}_{12}\left(x,h_4\right)\mathrm{d}x& h_4^2{\hat{f}}_2\left(h_4\right)& h_4\int x{\hat{f}}_{12}\left(x,h_4\right)\mathrm{d}x\\ {\hat{f}}_2\left(h_5\right)& \int x{\hat{f}}_{12}\left(x,h_5\right)\mathrm{d}x& h_5{\hat{f}}_2\left(h_5\right)& \int x^2{\hat{f}}_{12}\left(x,h_5\right)\mathrm{d}x& h_5^2{\hat{f}}_2\left(h_5\right)& h_5\int x{\hat{f}}_{12}\left(x,h_5\right)\mathrm{d}x\\ 1& {\bar{X}}_{n_{12}}& {\bar{H}}_{n_{12}}& {\displaystyle \frac{1}{n_{12}}}\sum _{i=1}^{n_{12}} X_i^2& {\displaystyle \frac{1}{n_{12}}}\sum _{i=1}^{n_{12}} H_i^2& {\displaystyle \frac{1}{n_{12}}}\sum _{i=1}^{n_{12}} X_i{H}_i\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ a_4\\ a_5\end{array}\right)=\left(\begin{array}{c}\int y{\hat{f}}_{02}\left(y,h_1\right)\mathrm{d}y\\ \int y{\hat{f}}_{02}\left(y,h_2\right)\mathrm{d}y\\ \int y{\hat{f}}_{02}\left(y,h_3\right)\mathrm{d}y\\ \int y{\hat{f}}_{02}\left(y,h_4\right)\mathrm{d}y\\ \int y{\hat{f}}_{02}\left(y,h_5\right)\mathrm{d}y\\ {\bar{Y}}_{n_{02}}\end{array}\right)$$

(68–73)

2.1.5. Performance Evaluation of the Models with Residual Analysis

Recall Equation (15):

$$\int m_{012}(x,h){\hat{f}}_{12}(x,h)\mathrm{d}x=\int y{\hat{f}}_{02}(y,h)\mathrm{d}y$$

(15)

Taking $h=H_j$, for $j\in \left\{1,\dots ,n_{02}\right\}$, and dividing by ${\hat{f}}_2\left(H_j\right)$, we get

$$\widehat{\stackrel{ˆ}{m}}\left(H_j\right)=\stackrel{\mathbb{ˆ}}{\mathbb{E}}\left(m_{012}(X,H)\mid H=H_j\right)$$

(74)

$$={\displaystyle \frac{\int {\hat{m}}_{012}\left(x,H_j\right){\hat{f}}_{12}\left(x,H_j\right)\mathrm{d}x}{{\hat{f}}_2\left(H_j\right)}}$$

(75)

$$={\displaystyle \frac{\frac{1}{n_{12}{b}_n}\sum _{i=1}^{n_{12}} {\hat{m}}_{012}\left(X_i,H_j\right)K\left(\frac{H_j-H_i}{b_n}\right)}{\frac{1}{n_2{b}_n}\sum _{i=1}^{n_2} K\left(\frac{H_j-H_i}{b_n}\right)}}$$

(76)

$$={\displaystyle \frac{n_2\sum _{i=1}^{n_{12}} {\hat{m}}_{012}\left(X_i,H_j\right)K\left(\frac{H_j-H_i}{b_n}\right)}{n_{12}\sum _{i=1}^{n_2} K\left(\frac{H_j-H_i}{b_n}\right)}}$$

(77)

as predicted values. Note that $\widehat{\stackrel{ˆ}{m}}\left(H_j\right)$ should be close to

$$\hat{m}\left(H_j\right)={\displaystyle \frac{\int y{\hat{f}}_{02}\left(y,H_j\right)\mathrm{d}y}{{\hat{f}}_2\left(H_j\right)}}$$

(78)

$$={\displaystyle \frac{\frac{1}{n_{02}{b}_n}\sum _{i=1}^{n_{02}} Y_{i}K\left(\frac{H_j-H_i}{b_n}\right)}{\frac{1}{n_2{b}_n}\sum _{i=1}^{n_2} K\left(\frac{H_j-H_i}{b_n}\right)}}$$

(79)

$$={\displaystyle \frac{n_2\sum _{i=1}^{n_{02}} Y_{i}K\left(\frac{H_j-H_i}{b_n}\right)}{n_{02}\sum _{i=1}^{n_2} K\left(\frac{H_j-H_i}{b_n}\right)}}$$

(80)

because of Equation (15).

The residuals can then be computed as

$$r_j=Y_j-\widehat{\stackrel{ˆ}{m}}\left(H_j\right), \mathrm{for} j\in \left\{1,\dots ,n_{02}\right\}.$$

(81)

and Mean Squared Error as

$$\mathbf{M}\mathbf{S}\mathbf{E}={\displaystyle \frac{1}{n_{02}}}\sum _{j=1}^{n_{02}} r_j^2={\displaystyle \frac{1}{n_{02}}}\sum _{j=1}^{n_{02}} {\left(Y_j-\widehat{\stackrel{ˆ}{m}}\left(H_j\right)\right)}^2$$

(82)

2.2. Demonstration Data Set

2.2.1. Destructive Sampling and External Tree Characterization

More than 40 candidate trees were initially identified through random selection on the LiDAR map of two unmanaged sites within the mixed wood stands of the Petawawa Research Forest (PRF) in Ontario, Canada (Figure 1). Of these, 18 trees were excluded due to decay or major structural defects, including forking, crooks, rotten branches, dead tops, and stem lean greater than 5°. The final sample consisted of 14 white pine aged (113 - 131 years) and 8 red pine trees (aged 70 – 128 years), which were manually felled using two approaches: (1) full-tree felling followed by delimbing on the ground (Figure 2a), and (2) delimbing the standing tree, lowering the branches to the ground with a crane, and subsequently felling the bole (Figure 2b). While the second method provided greater control and minimized damage, the first method resulted in breakage of approximately 5–10% of the branches.

After felling, all branches were organized to ensure that each limb remained intact and correctly oriented, and the bole and branches were positioned on a level surface for measurement but close to the stump. Stump characteristics, including stump height at the lowest and highest points and stump diameter were recorded. A reference line was then marked along the stem from the cut base to the tip (Figure 2c), and total tree length of the felled part was measured. Stump height was noted as the average height from the ground at the highest and lowest elevation. Total tree height is defined as the total of stump height and total tree length. The height to the base of the live crown (crown base height, CBH) is defined as the position of the lowest living whorl containing at least one live branch and with continuous whorls above it. Stem diameter was recorded at 1-m intervals up to the height of the first living branch.

For each whorl, starting from the last living branch to the tree apex, its position from the base of the stem was recorded. A sample of how the stem diameter at each recorded whorl was measured is shown in Figure 2d. Sum of the length between the branches of the whorl exactly 5 cm below the whorl was calculated as the diameter of the stem at the whorl. For every branch within each whorl, measurements included position relative to the reference line, branch length, branch diameter, branch inclination, and vitality status (dead or live). Finally, the stem was divided into 5-m sections for further bucking and evaluation. A small notch was made along the reference line at each 5-m mark (needed to keep track of the original orientation of each log within each tree for stem reconstruction, after which the stem was cut into logs and each segment was numbered.

To facilitate visualization and validation of the detailed measurements collected from the felled trees, a digital tree-reconstruction tool was developed in-house using OpenGL (Figure 3d,e). The tool visualizes the tree segments to the scale. Each tree was validated against the visualization model and bi-temporal lidar data that was flown before and after the harvest. Due either to the 50-cm width constraint during log CT scanning or to branch losses of 5 - 10% resulting from the felling process, four white pine and two red pine trees were excluded from further analysis.

Finally, a total of 58 white pine and 20 red pine 5m logs were sent for log reconstruction and knot characterisation. Because the height of the trees was variable, the number of logs and their lengths varied between trees. Over 97% of the logs had a length in the 4.5–5.5 m range and 2–7 logs were cut per tree (Table 1).

2.2.2. Knot Characterization

Following harvesting at the PRF, the segmented logs were transported to the Institut national de la recherche scientifique (INRS) in Quebec City for internal characterization. The knot measurement protocol was executed in two primary stages: initial X-ray Computed Tomography (CT) scanning, followed by automated feature extraction and measurement using the CT2Opti software suite.

X- Ray CT Scanning and Image Acquisition

Log sections were scanned using a Siemens Somatom Definition 128 AS+ medical CT scanner (Siemens-Healthineers, Erlangen, Germany) available at the INRS. To accommodate the scanner’s 2.1 m length limit, each 5 m log was halved (Figure 4a) and scanned in two passes (Figure 4b). CT images were reconstructed at 1 mm intervals along the longitudinal axis with a 2 mm slice thickness. This 1 mm overlap was intentional, to ensure maximum knot detection across consecutive frames [38]. Scanning parameters included an X-ray voltage of 120 kV, a current of 275 mAs, a pitch of 0.5, and a variable field of view yielding a 0.6 mm x 0.6 mm pixel resolution. Image reconstruction utilized the I70f Safir 3 filter, producing approximately 2,500 images per 2.5 m section.

The resulting Hounsfield units (gray tones) were utilized as a proxy for wood density [39]. While knots typically appear brighter due to higher density relative to clear wood [40], sapwood and wet pockets can exhibit similar density profiles, potentially obscuring contrast. This challenge was addressed using the specialized extraction pipeline detailed below.

Knot Extraction and Measurements with CT2Opti

Internal features were extracted using CT2Opti and Optitek 10, a specialized image-processing suite (please consult [38,41] for details) that combines morphological filtering, thresholding, and edge detection. The extraction followed a refined pipeline that distinguished real knots from false positives by requiring objects to originate near the pith and extend radially toward the log surface. The three-step pipeline can be described as:

• Candidate Identification: Where potential knot sections were identified in each 2D image
• 3D Reconstruction: Where individual sections were grouped into 3D objects by connecting adjacent pixels with similar characteristics across the image stack and
• Biological Filtering: Real knots were distinguished from false positives (e.g., wet pockets) by applying a knot model requiring objects to originate near the pith and extend radially toward the log surface.

This process enabled the detailed 3D reconstruction of diverse log profiles, ranging from large-diameter white pine butt logs with significant flare (Figure 4d) and low-tapered intermediary sections (Figure 4e) to complex crown logs with smaller diameters and higher knot densities (Figure 4f–g). High-resolution visualizations further demonstrate the characteristic anatomical expansion of knot size as branches grow from the pith toward the outer bark (Figure 4h).

Following manual review and digital reassembly along the pith-line, to reconstruct the complete tree architecture, high resolution metrics were measured for each sample with care. Metrics extracted by CT2Opti includes Total number of knots; Total knot volume (cm³) along with Average knot volume, Minimum knot volume and Maximum knot volume; Average knot area (cm²) along with Minimum knot area and Maximum knot area (found at the surface of the log), and total log volume (cm³). Table 2 summarises the measurements extracted from the felled white and red pine logs. In addition to these metrics, K/T index was also computed. K/T is defined as the proportion of knot volume over the total volume of the reconstructed tree, which provides insights of between tree or species variations [42].

2.3. Application of the Model – Estimating Knot Characteristics of White and Red Pines

2.3.1. Regression Strategy and Data Preparation

As the covariates varied between the lower and upper portions of the tree, the Regression of Misaligned Covariates (RMC) was estimated independently for the Clear Stem (CS) and the Living Crown (LC) across both species, white pine (WP) and red pine (RP). In this framework, trees were grouped by species. To evaluate model performance, a leave-one-tree-out cross-validation approach was adopted: for each species, one tree was sequestered as the evaluation (test) dataset, while the remaining pooled trees served as the training dataset.

To account for varying total heights among the pooled samples, the vertical coordinate was transformed into a relative altitude (RelAltitude) variable. This dimensionless index ranges from 0 to 1, where:

• For the CS: 0 represents the base of the bole and 1 represents the highest point of the clear stem.
• For the LC: 0 represents the transition from the clear stem and 1 represents the topping height.

Prior to modeling, every tree was validated by ensuring the total bole length matched the cumulative lengths of the CT scanned logs following the knot characterization process.

Variable Selection: The internal characteristics serving as response variables (Y) included Knot Volume (KnotVolume), Surface Knot Diameter (SurfaceKnotDiameter), and Knot Size (KnotSize)—the latter derived from the surface diameter assuming a circular knot geometry. The external covariates utilized for the LC included Log Diameter (LogDiameter), Branch Length (BranchLength), and Branch Diameter (BranchDiameter), while the CS model relied primarily on LogDiameter.

2.3.2. Model Application to Different Pine Species

As previously described, given a small set of felled trees, a Leave-One-Tree-Out cross-validation approach was employed to train the RMC model. Data partitioning was strictly controlled at the tree level rather than the individual knot level. This was to ensure independence between training and validation sets. In this demonstration, knot volume was chosen as the response as it connects statistical estimation directly to cumulative defect risk along the bole, allowing model performance to be evaluated in terms that are immediately relevant to log sorting and wood quality assessment. In addition, covariates that are directly estimable from LiDAR point clouds are chosen, viz., BranchLength and BranchDiameter within the live crown, and LogDiameter for the clear stem, to facilitate seamless integration of the proposed RMC framework with operational remote-sensing workflows and enable stand-level prediction of knot volume without reliance on destructive measurements.

Discrete knot volume predictions were integrated along the length of the bole to generate cumulative response profiles characterizing the rate of defect accumulation from the stump to the tree top. These profiles provide a spatially explicit basis for identifying optimal bucking points and log-quality thresholds.

2.3.3. Model Formulations and Software

Four distinct variations of the generalized model were implemented to evaluate different structural assumptions regarding the relationship between the response (Y), the covariates (X), and the relative height (H):

• norm: A parametric approach assuming (Y, X, H) follows a trivariate normal distribution.
• linear: A model assuming the regression function m(x, h) is linear in both x and h.
• quad: A model assuming m(x, h) follows a quadratic polynomial relationship in x and h.
• cubic: A model assuming m(x, h) follows a cubic polynomial relationship in x and h.

It is to be noted that despite the implementation of a cubic model, this study focuses only on norm, linear and quad models to maintain the intended scope.

Performance of the models was evaluated using the mean squared error derived in the section above. All modeling and statistical analyses were performed using R (v.2026.01.0; [43]).

3. Results and Discussion

3.1. Summary of the Demonstration Dataset

The demonstration dataset comprised mature timber, with most sampled trees exceeding 111 years of age; specifically, WP ranged from 111 to 131 years and RP) from 70 to 128 years, with only a single specimen (RP4) representing a younger age class at 70 years (Table 1). All samples adhered to a 50 cm diameter restriction necessitated by the physical aperture limits of the CT scanner. The WP specimens exhibited an average diameter at breast height (DBH) of 43.2 cm and an average total height of 27.54 m, while RP specimens averaged 40.0 cm in DBH and 25.22 m in height. Anatomical complexity varied significantly between the species: WP trees averaged 181 whorls per tree compared to 129 for RP. Out of a total of 2,104 recorded branches (1,590 WP; 514 RP), approximately 60% of them were identified as dead.

The comprehensive internal wood characterisation spanned a total bole length of 155 m for WP and 87 m for RP, resulting in the characterization of 3,316 and 1,150 knots, respectively (Table 2). Of these knots, over 90% of the surface diameters were below 8 cm, and 90% of knot sizes were below 20 cm² for both species. However, knot volume exhibited greater inter-species variation: 90% of the knots in RP were below 150 cm³, whereas the 90th percentile for WP extended up to 400 cm³ (Table 2).

3.2. Process of RMC to the Demonstration Dataset

Figure 5 illustrates the implementation of the RMC framework in R, predicting internal knot volume from externally measured log diameter across the entire bole of a sample RP specimen. It documents a two-step estimation strategy: estimate the one-covariate conditional means $m_{\mathrm{0,2}}\left(h\right)=E[Y\mid H=h]$and $m_{\mathrm{1,2}}\left(h\right)=E[X\mid H=h]$from training data, then solve the integral equations that recover the two-covariate regression $m_{\mathrm{0,1},2}(x,h)=E[Y\mid X=x,H=h]$under a chosen finite-dimensional model (normal, linear, quadratic, cubic, etc.). In this example, the training phase utilized three RP trees, comprising a total of 1,107 RelAltitude, 361 LogDiameter and 726 KnotVolume observations. LogDiameter (x) and KnotVolume (y), are not aligned in space, yet both exhibit distinct correlations with the shared spatial index, RelAltitude (h).

The resulting predicted surfaces demonstrate the framework’s ability to reconcile these disjointed variables into a unified 3D modeling space. Aside from the data range, the distributions of the training and test variables are largely consistent. While KnotVolume exhibits significant noise in both datasets, most observations remain below 70 cm³. Within this range, all candidate models, apart from the quad, successfully track the central density of the measured data across the tree’s relative height. Quantitative prediction performance, as indicated by the MSE, shows that the norm) and linear models yield results closely aligned with the measured data and the optimized m-hat estimator.

3.3. Application of RMC to Two Different Species

One sample tree from each of the species was randomly selected (viz., WP11 and RP5) to evaluate application of RMC to the two pine species. Both are matured timber of similar ages, of 125 (WP11) and 128 years (RP5), but exhibit distinct structural differences that test the model’s flexibility (Table 3). While RP5 reached a greater height (29.92 m) compared to WP11 (26.46 m), WP11 has a longer clear stem (10.3 m vs 7.75 m). In the living crown part, RP5 has higher number of branches as well as lower proportion of dead branches. As a result, RP5 has a higher number of knots, and almost twice the amount of knot volume compared to WP11. The higher variance in knot volume (up to 428.60 cm³) and the distinct K/T ratios between the two trees provide a comprehensive dataset to evaluate the predictive accuracy of the RMC framework across different pine species.

Table 4 details the partitioning used across the two primary physiological segments (clear stem and live crown). In the clear stem, the dataset maintains a high concentration of knot observations relative to the primary covariate (LogDiameter) to effectively model the spatial misalignment between disjointed measurements. Reflecting the greater morphological complexity of the upper bole, the living crown had significantly larger training sets, for instance, training data comprised up to 2,130 knot observations for white pine, and 620 for red pine, accounting for the increased density and volumetric variance of the branch base.

The resulting estimated regressions after applying the RMC framework, which characterize these spatial relationships for both species, are presented below and illustrated in Figure 6.

Estimated regressions:

(a)

White pine – Clear stem with LogDiameter:

$\hat{m}\left(x,h\right)={\hat{m}}_{012}\left(x,h\right)=\stackrel{\mathbb{ˆ}}{\mathbb{E}}\left(Y\right|X=x,H=h)=5.69\mathrm{*}{10}^{+01}+4.47\mathrm{*}{10}^{-02}x+2.59h (nrmal$)

${\hat{m}}_{012}\left(x,h\right)=7.82\mathrm{*}{10}^{+01}-5.93\mathrm{*}{10}^{-01}x+8.41h \left(linear\right)$${\hat{m}}_{012}\left(x,h\right)=1.68\mathrm{*}{10}^{+02}+1.49\mathrm{*}{10}^{+01}x-2.01\mathrm{*}{10}^{+02}h-3.79\mathrm{*}{10}^{-01}{x}^2+2.23\mathrm{*}{10}^{+02}{h}^2-3.51xh\left(quad\right)$

• (b) White pine – live crown with BranchLength

$\hat{m}\left(x,h\right)={\hat{m}}_{012}\left(x,h\right)=\stackrel{\mathbb{ˆ}}{\mathbb{E}}\left(Y\right|X=x,H=h)=6.53\mathrm{*}{10}^{+01}-5.28\mathrm{*}{10}^{-01}x-4.05\mathrm{*}{10}^{+01}h \left(normal\right)$${\hat{m}}_{012}\left(x,h\right)=-3.20\mathrm{*}{10}^{+01}+2.52x+3.41h \left(linear\right)$${\hat{m}}_{012}\left(x,h\right)=7.14\mathrm{*}{10}^{+01}-2.20x+1.59\mathrm{*}{10}^{+01}h+4.01\mathrm{*}{10}^{-02}{x}^2-3.76h^2-1.76xh \left(quad\right)$

• (c) White pine – live crown with BranchDiameter

$\hat{m}(x,h)={\hat{m}}_{012}(x,h)=\stackrel{\mathbb{ˆ}}{\mathbb{E}}(Y\mid X=x,H=h)=1.11\mathrm{*}{10}^{+02}-1.01\mathrm{*}{10}^{+01}$ (normal)

${\hat{m}}_{012}(x,h)=8.19\mathrm{*}{10}^{+01}-6.04x-2.85h$ (linear)

${\hat{m}}_{012}(x,h)=4.79\mathrm{*}{10}^{+01}+6.04\mathrm{*}{10}^{+01}x+4.60\mathrm{*}{10}^{+01}h-1.04\mathrm{*}{10}^{+01}{x}^2-2.25\mathrm{*}{10}^{+01}{h}^2-2.36\mathrm{*}{10}^{+01}xh$ (quad)

• (d) Red pine – Clear stem:

$\hat{m}\left(x,h\right)=1.96\mathrm{*}{10}^{+01}+4.83\mathrm{*}{10}^{-01}x+2.16\mathrm{*}{10}^{+01}h \left(norm\right)$${\hat{m}}_{012}\left(x,h\right)=-1.64\mathrm{*}{10}^{+01}+2.25x-3.71h \left(linear\right)$${\hat{m}}_{012}\left(x,h\right)=3.15\mathrm{*}{10}^{+01}-2.95x+8.87h+9.04\mathrm{*}{10}^{-02}{x}^2-3.86\mathrm{*}{10}^{+01}{h}^2+1.98xh \left(quad\right)$

• (e) Red pine – live crown with BranchLength

$\hat{m}\left(x,h\right)={\hat{m}}_{012}\left(x,h\right)=\stackrel{\mathbb{ˆ}}{\mathbb{E}}\left(Y\right|X=x,H=h)=-9.34+3.49\mathrm{*}{10}^{-01}x-2.74\mathrm{*}{10}^{+01}h \left(normal\right)$${\hat{m}}_{012}\left(x,h\right)=7.88\mathrm{*}{10}^{+01}-9.02\mathrm{*}{10}^{-02}x-3.18h \left(linear\right)$${\hat{m}}_{012}\left(x,h\right)=8.65\mathrm{*}{10}^{+01}-9.58\mathrm{*}{10}^{-01}x+1.19\mathrm{*}{10}^{+01}h+2.87\mathrm{*}{10}^{-03}{x}^2-1.24h^2+2.25\mathrm{*}{10}^{-02}xh \left(quad\right)$

• (f) Red pine – live crown with BranchDiameter

$\hat{m}\left(x,h\right)={\hat{m}}_{012}\left(x,h\right)=\stackrel{\mathbb{ˆ}}{\mathbb{E}}\left(Y\right|X=x,H=h)=2.11\mathrm{*}{10}^{+02}-3.86\mathrm{*}{10}^{+01}x-3.37\mathrm{*}{10}^{+01}h \left(normal\right)$${\hat{m}}_{012}\left(x,h\right)=1.08\mathrm{*}{10}^{+02}-2.40\mathrm{*}{10}^{+01}x+4.31\mathrm{*}{10}^{+01}h \left(linear\right)$${\hat{m}}_{012}\left(x,h\right)=8.17\mathrm{*}{10}^{+01}-2.03\mathrm{*}{10}^{+01}x-4.21\mathrm{*}{10}^{+01}h+1.74x^2+5.98\mathrm{*}{10}^{+01}{h}^2-2.52xh\left(quad\right)$

It should be noted that the m-hat estimator (Nadaraya-Watson estimator) provides a theoretical ceiling (a theoretical maxima and a benchmark) for predictive accuracy by utilizing known response (y) values. It represents the lowest possible error achievable for a given data structure. In practice, implementation for unknown timber relies on the proposed parametric models that utilize only the covariates (here they are externally measurable variables from LiDAR like log diameter, branch length and height).

In this application, m-hat performs well across all variable combinations, species and physiological sections justifying a good benchmark (Figure 6 and Table 5). It yields low or closest MSE to measured deviation (sd), for example 37.98 vs 38.71 and 37.95 vs 38.05 in the clear stem of WP11 and RP5 respectively, with conservative spread avoiding extreme overpredictions at high values. Similarly in the living crown, the model performance is comparable for both sets of covariates, BranchLength (38.37 vs 41.09 and 82.3 vs 83.92 for WP11 and RP5 respectively) and BranchDiameter (82.3 vs 83.92 and 37.96 vs 38.03 for WP11 and RP5 respectively), indicating its robustness in complex structural heterogeneity well. This also emphasizes its role as a strong operational benchmark for LiDAR-based standing-tree wood-quality estimation.

Within the clear stem section, the predictive performance of the proposed models varied between the species (Table 5 and Table 6). Both the linear and normal models closely approximated the trajectory of the m-hat benchmark across the bole, with linear model having a slightly higher bias at higher altitudes in WP11 compared to RP5 (Figure 6). Comparing their predictive performance, both models have lower bias and consistently near-parity ratio against that of m-hat, establishing their reliability for predicting knot volume from external LogDiameter measurements. Conversely, the quadratic (Quad) model exhibited a consistent upward bias at higher RelAltitude and five times higher error compared to m-hat. This sensitivity to extreme values was particularly pronounced in the morphologically complex white pine (WP11) relative to the more uniform red pine (RP5).

Predictions in the crown section show higher overall error across all models due to increased data variance, yet the hierarchy of performance remains consistent (Figure 6, Table 5 and Table 6). For WP11,BranchLength as well as BranchDiameter all models have near-parity with the m-hat, especially normal performed within 10% of m-hat, but showed higher mean bias in estimating KnotVolume when using BranchLength. For RP5, the quadratic model performed surprisingly well in the crown using BranchDiameter (41.92 vs. 37.96 for m-hat and low mean bias), suggesting that for this specific species and section, a slightly curved relationship better captures the volume of larger branch bases. However, linear model performs well with BranchLength as the predictor closest to m-hat with a comparatively low mean bias (Table 5 and Table 6).

3.4. Management Implications – Turning Data into Decisions for Bucking and Log Sort

Profiles of the cumulative estimated knot volume along the bole height for the proposed and benchmark models (WP11 and RP5) are illustrated in Figure 7. Vertical markers indicate the heights at which 25%, 50%, and 75% of the total knot volume are exceeded. These profiles provide a high-resolution map of internal wood quality for individual standing trees, identifying transition points and accumulation rates that allow foresters to move beyond external visual assessments toward more precise, volume-based grading.

For both species, the 25% threshold identifies the portion of the bole most likely to yield high-quality butt logs with minimal knot impact. When this threshold lies above the clear stem height, the first log can be confidently assigned to sawlog or veneer classes. The 50% threshold marks a transition zone where knot influence becomes substantial; logs cut above this height are increasingly likely to require downgrading to lower-grade sawlogs and lumber. Finally, the 75% threshold suggests that the upper bole sections contribute little additional high-quality material and are better directed to pulp or biomass streams.

Except for the quadratic model in WP11, all proposed models effectively captured the accumulation trends observed in the benchmark and measured data profiles. Although the clear stem of WP11 (10.3 m) is longer than that of RP5 (7.75 m), WP11 accumulates the first 25% of its knot volume at a significantly faster rate. Beyond this point, accumulation in WP11 becomes gradual, whereas RP5—with its more complex branch structure—adds knot volume rapidly after the 50% threshold. Specifically, while WP11 reaches the 25% mark at 10.7 m (near the clear stem top), RP5 does not reach this same threshold until 13.8 m.

Species-specific differences in these threshold positions provide a framework for differentiated sorting strategies. WP11’s later accumulation of most of its knot volume offers a greater proportion of high-grade lower logs. Conversely, the earlier accumulation trends in RP5 imply that a more conservative bucking height is required for premium products. Ultimately, these cumulative curves allow managers to predefine bucking heights that balance value recovery and defect risk, improving consistency and efficiency in industrial log sorting.

4. Discussion

This study introduces a generalized statistical framework, Regression of Misaligned Covariates, for reconciling covariates and responses that are measured at different spatial locations - a long-standing challenge in wood-quality modelling and in the integration of multi-source forest data. By combining the Law of Total Expectation with kernel-based estimation of marginal relationships, the RMC approach provides a principled way to recover conditional mean structure even when the underlying datasets are spatially misaligned. This directly addresses the “change of support” problem that arises when LiDAR-derived external attributes and destructively sampled internal knot properties are not observed at the same coordinates.

Applying RMC to eastern white pine and red pine demonstrates the flexibility of the framework across contrasting physiological architectures. Red pine, with its uniform crown structure, provided a stable baseline, while eastern white pine with its irregular branching and susceptibility to structural disturbances, served as a rigorous test of robustness. Across both species, RMC produced smooth, biologically interpretable patterns of knot volume along the bole, capturing species-specific differences in the rate and height at which knot volume accumulates.

The study also demonstrates simplicity of its application and operational practicality. The results highlight that the method does not require species-specific parameterization or complex covariate sets, history or stand conditions. In fact, the model performed well using only height and a single LiDAR-derivable covariate. Notably, the framework produced stable estimates even when trained on a very small dataset, for instance only four red pine trees were used to train and predict RP5 knot characteristics, suggesting that RMC can extract meaningful structure from limited data, though performance will naturally improve as larger datasets become available.

The cumulative knot-volume profiles derived from RMC also illustrate the direct management value of the framework. Expressing knot accumulation as a function of height provides a clear link to operational decisions such as bucking, log sorting, and product allocation. Thresholds marking where 25%, 50%, and 75% of total knot volume occur allow managers to identify the heights at which log quality begins to decline, assess the likely grade of the first and second logs, and anticipate the proportion of the stem suitable for sawlog, dimension lumber and/or studwood, or pulp. Because these thresholds are expressed as cumulative proportions, they can be compared across species, stands, and management regimes, making the approach broadly applicable in both industrial and silvicultural contexts.

Beyond wood-quality modelling, RMC represents a general framework for any application involving misaligned covariates and responses. Many ecological, environmental, and remote-sensing datasets share this structure—for example, soil properties measured at sparse points paired with canopy metrics from LiDAR, or wildlife observations paired with environmental covariates collected at different scales. The RMC formulation is therefore not limited to forestry; it provides a transferable solution for a broad class of misalignment problems where traditional regression approaches fail.

There are, however, several avenues for further work. First, the choice of kernel bandwidth $\left(b_n\right)$plays a central role in shaping the smoothness and accuracy of the marginal estimates that feed into the RMC solution. A careful, data-driven selection of $b_n$—potentially through cross-validation or adaptive bandwidth strategies—will be important for optimizing performance in different forest conditions. Second, the present study demonstrates the framework using two covariates (height and one structural attribute) to illustrate the core idea. Extending RMC to three or more covariates $\left(H,X_1,\dots ,X_p\right)$is a natural next step and would allow the integration of richer LiDAR-derived metrics such as branch angle, crown width, or local curvature of the stem. Finally, although the method performed well on small datasets, broader validation across larger samples, additional species, and diverse stand structures will help refine its operational reliability.

Overall, RMC provides a portable and adaptable modelling framework that can be applied to any species or region without requiring species-specific calibration. Because the covariates used here are directly measurable from terrestrial or airborne LiDAR, the approach enables non-destructive prediction of internal wood quality in standing trees. This positions RMC as a practical tool for pre-harvest planning, value-based harvesting, and long-term monitoring of stand development - supporting more precise and sustainable forest management.

5. Conclusions

This study presents RMC as a general and scalable solution to the persistent problem of misaligned covariates and responses in forestry and beyond. By validating the framework across two physiologically contrasting conifer species, we demonstrate that reliable, biologically interpretable estimates of knot volume can be recovered even from small datasets and using only minimal LiDAR-derived covariates. The resulting cumulative profiles translate directly into operational guidance for bucking, sorting, and value-based harvesting, while also enabling non-destructive assessment of standing-tree wood quality. Because RMC is not tied to species-specific parameters and can incorporate additional covariates as needed, it is readily portable across regions, forest types, and remote-sensing platforms. As multi-source data become increasingly central to forest monitoring, RMC offers a flexible pathway for integrating external structural measurements with internal wood-quality attributes, supporting both pre-harvest decision-making and long-term sustainable management.