**"Tomographic image reconstruction using an interpolation method for tree decay detection,"**

*BioRes.*9(2), 3248-3263.

#### Abstract

Stress wave velocity has been traditionally regarded as an indicator of the extent of damage inside wood. This paper aimed to detect internal decay of urban trees through reconstructing tomographic image of the cross section of a tree trunk. A grid model covering the cross section area of a tree trunk was defined with some assumptions. Stress wave data were processed beforehand to obtain the propagation velocity and the coordinate values. An image reconstruction algorithm for detecting internal decay was proposed based on an interpolation method, which estimated the velocity values of unknown grid points by utilizing the values of the surrounding points. To test the effectiveness of this method, *Cinnamomum camphora* tree samples were selected and tested using a stress wave tool. The area, positions, and extent of decay in the representative samples were displayed in tomographic images constructed by the interpolation method, and the results demonstrate the performance of the method.

Download PDF

#### Full Article

**Tomographic Image Reconstruction Using an Interpolation Method for Tree Decay Detection**

Hailin Feng,^{a,}* Guanghui Li,^{a} Sheng Fu,^{a} and Xiping Wang ^{b}

Stress wave velocity has been traditionally regarded as an indicator of the extent of damage inside wood. This paper aimed to detect internal decay of urban trees through reconstructing tomographic image of the cross section of a tree trunk. A grid model covering the cross section area of a tree trunk was defined with some assumptions. Stress wave data were processed beforehand to obtain the propagation velocity and the coordinate values. An image reconstruction algorithm for detecting internal decay was proposed based on an interpolation method, which estimated the velocity values of unknown grid points by utilizing the values of the surrounding points. To test the effectiveness of this method, *Cinnamomum camphora* tree samples were selected and tested using a stress wave tool. The area, positions, and extent of decay in the representative samples were displayed in tomographic images constructed by the interpolation method, and the results demonstrate the performance of the method.

*Keywords: Stress wave; Tomography; Internal decay; Nondestructive testing; Trees*

*Contact information: a: Information Engineering School of Zhejiang Agriculture and Forestry University, Hangzhou, China; b: USDA Forest Service, Forest Products Laboratory, Madison, Wisconsin, USA;*

** Corresponding author: sealinfeng@gmail.com*

**INTRODUCTION**

Internal decay in trees reduces the value and quality of wood and causes economic loss to forests and raw wood materials that are used to manufacture various wood products. Detecting internal wood decay and assessing a tree’s safety using nondestructive methods has received increasing attention in recent years in both research and application.

The acoustic method has been found very effective for detecting internal decay of urban trees. Stress wave velocity has been traditionally regarded as an indicator of the presence and extent of internal wood decay (Tomikawa *et al.* 1986; Auty and Achim 2008; Roohnia *et al.* 2011). Generally, stress waves travel more slowly in decayed or deteriorated wood than in sound wood, as decayed wood could cause a prolonged wave transmission time in measurements (Wang 2013). When a tree suffers from decay or other internal damages, stress wave propagation velocity changes according to the change of its physical conditions. Therefore, when conducting the acoustic measurement, the internal condition of a tree can therefore be detected by analyzing the stress wave data.

In tree inspection practice, a single-path stress wave approach is commonly used to measure the wave transmission time through the cross section of a tree. The wave transmission time can then be used to make a preliminary diagnosis of whether the tree trunk is sound or not. This approach is useful for preliminary screening of urban trees and identifying trees that could potentially have decay or other internal defects. However, to make a better assessment of the true physical condition of a tree, more scientific information relating to the internal condition of trees needs to be acquired and analyzed, which leads to development of more sophisticated tomography tools with high accuracy (Wang *et al.* 2009).

Ultrasonic computed tomography was first used to inspect wooden poles by Tomikawa *et al. *(1986). Since then, acoustic tomography technology has been further investigated for its use in the field of tree inspection (Rabe *et al.* 2004; Gilbert and Smiley 2004; Socco *et al. *2004; Yanagida *et al.* 2007; Deflorio *et al.* 2008; Lin *et al.* 2008; Wang *et al.* 2009; Schubert *et al.* 2009; Brazee *et al.* 2011; Li *et al.*2012). The tomography technologies currently being used for tree inspection are mostly based on time-of-flight (TOF) of stress waves or electrical resistance of electromagnetic waves.

Different inversion algorithms have been tested in acoustic tomography. Filtered back-projections (FBP) can be used to perform an inverse transformation to obtain the slowness function from the TOF by assuming that stress waves propagate in straight lines (Brancheriau *et al.* 2012). The simultaneous iterative reconstruction technique (SIRT) (Lakshminarayanan and Lent 1979) is an iterative root mean square optimization (RMSO) (Socco *et al.* 2004) that calculates inverse solutions of a trunk model based on measured TOF data, which is essentially an algebraic reconstruction technique suffering from an ill-conditioned problem. Assuming that a transverse cross-section of wood is isotropic, Brancheriau *et al. *(2008) designed a reconstruction algorithm using a linear approximation of the forward problem. Nicolotti *et al.* (2003) evaluated three different types of tomography methods for calculating tomograms of decay trunks in urban trees. The reconstruction algorithm, the acoustic frequency, and the numbers of sensors all have influence on the resolution of images obtained by stress wave based acoustics (Divos and Divos 2005).

Although acoustic tomography has been found useful in detecting internal decay of urban trees, the potential of this technology has not yet been fully explored. Complicated factors still limit the effectiveness of acoustic tomography in field applications. The reconstructed tomographic images suffer from significant alteration due to the property of projections being a poor approximation of the object area. Large variations exist in stress wave propagation velocity, which may be caused by the physical conditions in different angles and the expected internal decay of a tree. Furthermore, anisotropy of wood is also a complicated factor, though an anisotropy travel time data correction factor was proposed for detection of decay in trees (Maurer *et al. *2006).

Based on the assumption that stress wave propagates as straight lines in the cross section of a tree, this paper specifically investigated the effectiveness of using an interpolation method to reconstruct a tomographic image from stress wave TOF data. A novel image construction method was proposed using a spatial interpolation technology for tree decay detection. Our main consideration is to predict the values of the unknown points by processing the surrounding variables with meaningful values within the same region. It was hypothesized that every point in the cross section of a tree trunk is related to other points and the adjacent points are more closely related than distant points. Consequently, the near points are usually assigned higher weight than distant points, which is the basic premise behind our image reconstruction process. With all the values of estimated nodes, the internal image of the cross section of a tree can be constructed using a certain color scale.

In this study, *Cinnamomum camphora* trees were tested to examine the effectiveness of our method. The tomographic images of the tested cross sections were reconstructed using the spatial interpolation method. The area, positions, and extent of the decay in the representative samples were displayed by the proposed method, and the results demonstrate the performance of our method.

**Theoretical Consideration**

In the cross sectional region of a tree trunk, first a discrete model is defined that represents the reconstruction area covering the cross section, and then the propagation velocity of a stress wave in the region is estimated. A grid model is then defined, and an interpolation scheme is used to estimate the values of all grid nodes. A spatial interpolation scheme can predict the values of the unknown points by processing the surrounding variables with meaningful values within the same region, which is often used to transfer the measured data on discrete points to a continuous data surface. It is supposed that the similar possibility is high if the unknown points are close proximity to the known points.

It is also considered that the TOF of a stress wave propagating across a tree trunk is related to the physical condition of the tree, and the propagation velocity has a spatial correlation within the region of interest. Determining the value of wave propagation velocity at every point in the cross section is the key in assessing the internal condition of the tree. The propagation velocity of a stress wave is treated as a variable distributed in the region, and the features of internal condition are reflected by the values of the variable.

A tree trunk under testing is considered as a cylindrical body and its cross section an ideal circle. The cross section of the trunk is defined as the target area Ω. The following analysis will be limited to this target area. The coordinates of points in a space can be represented as and the propagation velocity of a stress wave can be taken as the observed value of the point. The observed values of the points in area Ω vary due to the anisotropic nature and internal conditions of wood material. It was regarded here as a regionalized variable. There is an autocorrelation in some degree between variables at the point *x*with those at point (*x+h*). The autocorrelation depends on the distance between two points and the variable properties. The variable in area Ω should satisfy the following conditions:

(1) Space limitation. The variable is limited to the cross sectional area in interest, which is defined as a regionalized geometric domain.

(2) Continuity. It can be described by the variograms of the regionalized variable.

(3) Space correlation. The variables in certain region have a degree of space correlation, which will extend beyond the region.

**Methodology**

First, a discrete model was created using a grid to represent the cross-sectional area of a tree trunk of interest. It was assumed that stress waves propagate from a source sensor to receiving sensors and pass through the cross section in a straight path. The coordinates and propagation velocity of mesh points in the cross section can be expressed as , where and are the coordinates of the intersection points of the stress wave, and is the propagation velocity measured through different sensors. A grid model was built to define the position of all sensors. The distribution of the intersection points of stress waves in the Cartesian coordinate system in Fig. 1. represent sensors, and the total number of sensors is represented by *t*. Given that the number of every sensor is *n*, *t *and* n *are positive integers, and , are the arc lengths which are measured in the testing process. The radius of the cross-section can be expressed as *,* where *t* is a positive integer. Let *A _{1 }*be the reference point and the other sensors

*A*to

_{2}*A*be distributed on the circumference in a clockwise direction.

_{t}**Fig. 1.** Coordinate system of the receiving sensors.

Let be the clockwise angle between any sensor *A _{n}* and the reference sensor

*A*. Because it is given that and , where

_{1}*t*and

*n*are both positive integers, so when , one has:

(1)

wherein , . and when n=1.

The equation of the cross-section circles of the testing wood can be shown as

(2)

Accordingly, the parameter equation of sensor *A _{n}* with coordinates (

*X*) can be expressed as:

_{n}, Y_{n}(3)

Let *R _{n }*be the radius offset,

*L’*be the length offset, and both terms are used to adjust the positions of the sensors. In the conditions of and , the deviation angle can be computed as:

(4)

In the equation, the position of the sensors can be adjusted by *R _{n}* and

*L’*. The coordinates of any sensor

*A*can be shown as:

_{n}(5)

Next, the coordinates of all sensors are determined. The intersection points of wave propagation paths in a tree cross section now need to be determined. By drawing a line connecting two sensors specified by the coordinate pairs, the coordinates of all intersection points can be ascertained. If two stress wave lines intersect at one point inside the circle but not the end point (position of a sensor), then it is considered an intersection point. There are three possibilities to confirm whether two given stress wave lines have intersection points.

(a)

(b)

(c)

**Fig. 2.** (a) Standard intersections between the stress wave lines; (b) No intersections between the stress wave lines; (c) Non-standard intersection between the stress wave lines

Take the measurement with 8 stress wave sensors as an example. Figure 2(a) indicates possible intersection points for two stress wave lines among four sensors (*i.e.,* ) which is shown by the dotted lines. There is no intersection point of two stress wave lines among four sensors in Fig. 2(b) depicted as the two dotted lines. In Fig. 2(c), there is one intersection point at for the two stress wave lines among three sensors , which technically is not considered an intersection point. As a result, the coordinates of the sensors and the intersection of two lines are used as the input of our interpolation method and can compute the coordinates of all intersections.

Following the definition of the grid model, we suppose that *x* is any point in the target region. Here, the region is the cross section of a trunk, and *z*(*x _{i}*) is the measured value. The code number of measured value in target area is defined as

*N*. The estimated value

*z*at point

^{*}(x_{0})*x*can be estimated by the

_{0}*N*known values within the area, and the

*N*known values

*z*(

*x*)(

_{i}*i=1,2,…,N*) are measured from the

*N*variables in the region limited to the cross section of the tree trunk.

Using the interpolation method, linear combination of those known values is used to express the estimated value* z ^{*}(x_{0})* as in the following the formula (Beers and Kleijnen 2004),

(6)

Where is a weighting coefficient, which determines the influence of the known values* z(x _{i})* on the estimated value

*z*. Obviously, the ability to estimate the values depends on the calculation and selection of the weighting coefficient.

^{*}(x_{0})There are two conditions that need to be met when trying to find the weighting coefficients. The first condition is that *z ^{*}(x_{0}) *should be an unbiased estimator. The formal description of the condition can be expressed as:

(7)

By combining Eqs. 6 and 7, one obtains:

(8)

The second condition is that the estimator should be optimal, which means the sum of the square of the deviation between the estimated value *Z ^{*}(x_{0}) *and the actual value is minimum, which can be represented as follows.

(9)

Then one has the following equation.

(10)

Now the task transfers into an optimization problem to minimize the target function within the parameter constraint. The Lagrangian method can be used to solve the problem, and a function is constructed, where is the Lagrangian operator,

(11)

where, (*x _{i}*,

*x*) is the semi-variogram value and the parameter

_{j}*h*is the distance between grid nodes

*x*and

_{i }*x*. Combining the Eqs. 9 and 11, one gets the n+1 order linear equations.

_{j}(12)

The linear equations can be expressed by a matrix form:

(13)

where is the semi-variogram value (*x _{i}*,

*x*). The weight parameters can then be obtained by solving the simultaneous Eq.12 and matrix 13. With those weight parameters, the estimated values of the nodes in the cross section can be deduced as the sum of the products of the weight parameters and the values of the input nodes with Eq. 6. Then next estimated value was calculated and the steps above were repeated until all estimations were completed. With all the values of estimated nodes and the coordinates, the internal image of the cross section of a trunk can be constructed with a certain color scale.

_{j}**EXPERIMENTAL**

Four *Cinnamomum camphora* tree samples with different diameters were selected as specimens to test the feasibility of our method. The first test sample was a standing tree with a diameter of 79 cm, located on the campus of Zhejiang Agriculture and Forestry University in Huangzhou, China. Stress wave measurements were conducted on the tree at a height of 1.2 m above the ground. The second sample was a rotted tree stump also located on the campus, with a diameter of 50 cm at the height of 30 cm above the ground. The internal trunk decay of the tree stump can be visually seen which enabled a direct comparison with our tomography results. The third and fourth test samples were 50-cm long trunk sections cut from two felled* *trees, with a diameter of 36 and 44 cm respectively. Stress wave experiments on these two samples were conducted in the laboratory. The end cross sections of these two samples can be visually observed and the tomographic images can be directly compared with the images of the cross sections. Stress wave measurements were repeated multiple times on each sample to examine the repeatability of the measurement data.

A stress wave device named WoodPecker was developed at the Laboratory of Intelligent Measurement Technology of Zhejiang Agriculture and Forestry University. The device consists of 12 sensors, a signal processor box, and a hand-held steel hammer. During measurement, 12 sensors were mounted around the tree trunk (on-site) or trunk section (in laboratory) through nails and equally spaced around the circumference. Stress waves were generated by impacting individual sensors using the hammer. A complete measurement on each sample required sequentially tapping each sensor until a TOF data matrix was obtained and recorded. When a sensor is hit, an impulse energy is introduced into the tree trunk generating low frequency (less than 7 kHz mostly) and high energy stress waves. All the remaining sensors receive the signals. The signal processor analyzes the signals and then determines the TOF between two sensors. For each test sample, the circumference and distances between sensors were measured using a tape measure and a caliper. This information was used to map the approximate geometric form of the cross section and determine the coordinates of points in the grid model. Upon completion of stress wave measurement on a sample, a tomographic image was constructed for the cross section tested using the software which was developed based on the spatial interpolation method discussed in earlier section.

**RESULTS AND DISCUSSION**

In the tomographic images constructed for the test samples, five color schemes (dark blue, light blue, yellow, orange, and red) were used to represent the internal conditions of the cross section tested. Dark blue represents the high stress wave velocity, high density, and sound wood; light blue, yellow, and orange represent areas of relatively low stress wave velocity, low density, and unsound or potentially decayed wood, respectively; red represents the lowest stress wave velocity with severe decay or void. The color scale is expanded between the 100% and lowest velocity.

In principle, a stress wave tomogram shows the distribution of stress wave velocity in the cross section of a tree. The color presentation of the velocity zones demonstrates the differences in the ability of the wood to transmit stress waves. The ability of the wood to transmit stress wave signals within a cross-section strongly correlates with its modulus of elasticity and density.

**Fig. 3.** Stress wave velocities of different orientations in standing tree of sample 1

Figure 3 shows the velocity trends in different angles in the cross section in test sample 1 (standing tree). It can be seen from the data that the velocity changed as the wave propagation path changed across a cross section. Every velocity trend line in Fig. 3 implies the velocity from a certain source sensor. The orientation represents the different direction to the source sensor in the circumference from tangential direction to radial direction. It has been shown that in sound wood there exists a second order polynomial relationship between the velocity ratio and wave propagation path angle (Li 2013). Compared with that in sound wood, the velocity patterns for Fig. 3 have a different representation. The stress wave velocity fluctuated in different directions, which is in accordance with the velocity pattern in wood with slight or early decay. A precise diagnosis of a tree’s condition can be reached by reading the tomogram which relates the velocity value and velocity pattern to the internal condition.

**Fig. 4.** Spatial distribution of the stress wave velocities and the contour map in standing tree of sample 1

**Fig. 5.** The reconstruction tomogram and the testing photograph of the tree sample 1

The cross-section maps of sample 1 with measured stress wave velocity are shown in Figs. 4 and 5. The values of stress wave propagation velocity are shown in Fig 4a as a spatial distribution in the cross section, which are the foundation for tomographic image reconstruction. Based on the grid model developed, the coordinates of all the intersection points and the velocity are presented. Every point in the map is expressed as , where and are the coordinates of the intersection points of the stress wave paths, and is the propagation velocity at each point.

Taking the data in Fig. 4(a) as input, calculations were made to estimate the data of unknown points by solving the equations. After the estimated data and coordinates of the unknown points in the cross section were derived, a contour map was developed for the tree sample (Fig. 4b). The color scale of the contour map is shown on the right side. It should be pointed out that the color code used in our image reconstruction is similar with that of the other inversion method, although the particular colors are not identical.

It was then necessary to remove the zigzags on contours, which produced a smoother reconstruction map without losing the original structures of the contours. The reconstruction map is shown in Fig. 5(a) with the same color scheme, and the color scale of the tomogram is also shown on the right side of the map. It is implied from Fig. 5(a) that there was a slight decay in the center of the cross section of the tree, and the area and position of the decay is displayed in the tomogram. The extent of the decay in the sample is represented by different colors ranging from yellow to red, which is also an indicator of stress wave propagation velocity. From the velocity patterns in Fig. 3, it was speculated that the physical conditions were heterogeneous and there may have been slight decay present. Compared with the velocity patterns, the internal condition of wood was revealed by the stress wave maps, and the same conclusion can be drawn from our mapping results.