Research on Tree Pith Location in Radial Direction Based on Terrestrial Laser Scanning

Obtaining the direction of a diameter line through the tree pith is the basis of effective sampling by a micro-drill resistance instrument. In order to implement non-destructive tree pith location in the radial direction, the geometric property of tree pith, the longest chord through the tree pith on the cross-section will bisect outer contour circumference, as first proposed and proven in this paper. Based on this property, a non-destructive tree pith radial location method based on terrestrial laser scanning was developed. The experiments of pith radial location were made on the tree discs and the error of location is less than 1.5% for cross-section shape closed to ellipse on four tree species. The geometric property and location method of the tree pith in this research would play an important role in studying the growth process of standing trees, obtaining processed wood properties, and estimating tree age.


Introduction
The tree trunk consists of bark, cambium, xylem, and tree pith [1]. The pith is ground tissue formed by the ground meristem, whereas the primary xylem is formed by the procambium. As the foundation of forest survey and trunk analysis, pith location plays an important role in measuring the parameters of tree ring, studying the growth process of trees, obtaining wood properties, and estimating tree age [2][3][4].
In order to evaluate the actual state of the tree or wood structure, there are usually four pith-related ways. (1) The method based on image processing: Chen Jie et al. scanned the images of tree disk horizontally and vertically, then detected the tree piths one by one in horizontal and vertical coordinates [5]. (2) The method based on CT [6][7][8][9][10][11][12][13][14]: Alkan et al. got the CT images of the tree disks, then some image processing methods were adopted to locate the tree pith and reconstruct the 3D images [15]. However, one limitation of the two above-mentioned methods is the need to cut down trees. (3) With their method based on ultrasonic wave, Perlin et al. placed twenty fixed ultrasonic receivers at equal intervals along the tree disc, the pith is located at the convergence point of the maximum sum of all ultrasonic velocities [16]. Results showed that the distance error between the final result and the actual pith was 0.5 cm. (4) The method based on using a micro-drilling resistance instrument: Perlin et al. pointed out that this method could be widely used in tree ring parameter measurement and tree age estimation due to the relative lower cost [16].
However, one of the main problems of micro-drilling resistance instruments is that drill bits do not reach the pith due to difficult alignment when drilling or pith position deviating from the cross-section center [17,18]. When coring is off the pith, several holes must be implemented till obtaining satisfactory results [16]. The potential solution is the pith position in radial direction is estimated before drilling. The tree pith location in radial direction is to determine the direction of a diameter line, passing through the tree pith. The pith radial location outside of the cross-section, is more challenging when the cross-section 2 of 14 is visible, but this method would work for micro-drilling resistance instruments coring in standing tress and processed wood.
Nevertheless, related researches on radial tree pith positioning are relatively few. Zhang Tongwen et al. designed a tree pith location calibration instrument to ensure that the drilling position of micro-drill resistance instrument can pass through the pith. The working principle is that the pith is located in the bisector of cross-section outer contour circumference of trunk [19]. However, the result of this method is accurate only when the pith is located in the geometric center of the cross-section. When the pith is off the center of cross-section, the bisectors of the circumference obtained by this method may not pass through the tree pith.
The purpose of this paper is to present a pith radial location method outside of the cross-section [20]. First, a hypothesis about the geometric properties of tree pith on the cross-section was proposed for the first time, and then the hypothesis was verified. Second, in order to easily compare the location result of this method with the real pith position, tree discs were employed to obtain point cloud of cross-section by 2D LiDAR, then tree pith location in radial direction was realized after series of data processing.

Hypothesis
Cambium only grows one layer around the pith each year, named a ring. Therefore, the pith is usually located in the center of the trunk. We can conclude that the number of the rings, which the diameter lines pass through, must be the largest on the cross-section of the trunk. Combining the above-mentioned researches, the hypothesis about the geometric properties of tree pith is presented in this paper [21][22][23][24][25], that is, the longest chord passing through tree pith on the cross-section can bisect outer contour circumference.

Hypothesis Verification
In order to testify the hypothesis of the geometric properties of the tree pith, the digital images of tree discs were imposed into AutoCAD [26,27]. Then, auxiliary lines were drawn and measured. Finally, the statistical analysis method was used to verify whether the results met the hypothesis.

•
Image acquisition of tree discs The discs of four different species of trees were adopted as the research samples in this paper, including camphorwood, chestnut, China fir and pine. These tree species are widely distributed in China and are common wood raw material. Accordingly, a total of 12 transverse-cutting sample discs were employed for each species and the digital images were mostly obtained from the internet and a small part were obtained by the camera. Some samples are shown in Figure 1. The pith radial location outside of the cross-section, is more challenging when the crosssection is visible, but this method would work for micro-drilling resistance instruments coring in standing tress and processed wood. Nevertheless, related researches on radial tree pith positioning are relatively few. Zhang Tongwen et al. designed a tree pith location calibration instrument to ensure that the drilling position of micro-drill resistance instrument can pass through the pith. The working principle is that the pith is located in the bisector of cross-section outer contour circumference of trunk [19]. However, the result of this method is accurate only when the pith is located in the geometric center of the cross-section. When the pith is off the center of cross-section, the bisectors of the circumference obtained by this method may not pass through the tree pith.
The purpose of this paper is to present a pith radial location method outside of the cross-section [20]. First, a hypothesis about the geometric properties of tree pith on the cross-section was proposed for the first time, and then the hypothesis was verified. Second, in order to easily compare the location result of this method with the real pith position, tree discs were employed to obtain point cloud of cross-section by 2D LiDAR, then tree pith location in radial direction was realized after series of data processing.

Hypothesis
Cambium only grows one layer around the pith each year, named a ring. Therefore, the pith is usually located in the center of the trunk. We can conclude that the number of the rings, which the diameter lines pass through, must be the largest on the cross-section of the trunk. Combining the above-mentioned researches, the hypothesis about the geometric properties of tree pith is presented in this paper [21][22][23][24][25], that is, the longest chord passing through tree pith on the cross-section can bisect outer contour circumference.

Hypothesis Verification
In order to testify the hypothesis of the geometric properties of the tree pith, the digital images of tree discs were imposed into AutoCAD [26,27]. Then, auxiliary lines were drawn and measured. Finally, the statistical analysis method was used to verify whether the results met the hypothesis.  Image acquisition of tree discs The discs of four different species of trees were adopted as the research samples in this paper, including camphorwood, chestnut, China fir and pine. These tree species are widely distributed in China and are common wood raw material. Accordingly, a total of 12 transverse-cutting sample discs were employed for each species and the digital images were mostly obtained from the internet and a small part were obtained by the camera. Some samples are shown in Figure 1.

Camphor wood_1
Camphor wood_2 Camphor wood_3 Camphor wood_4 Camphor wood_5 Camphor wood_6 Chestnut_1 Chestnut_2 Chestnut_3 Chestnut_4 Chestnut_5 Chestnut_6   Sample images were shown in Figure 2. After sample images were imported into AutoCAD, the pith positions were marked manually. Then, as shown in Figure 3, many relatively long chords through the pith were assumed, measured and compared. Finally, the longest chord passing through the pith on the tree disc was obtained.

Perimeter measurement
In order to measure the length of every part of the outer contour of the tree disc, split by the longest chord through the pith, polylines were respectively made on the two parts, based on the principle of cyclotomy.   Figure 2. After sample images were imported into AutoCAD, the pith positions were marked manually. Then, as shown in Figure 3, many relatively long chords through the pith were assumed, measured and compared. Finally, the longest chord passing through the pith on the tree disc was obtained.  Sample images were shown in Figure 2. After sample images were imported into AutoCAD, the pith positions were marked manually. Then, as shown in Figure 3, many relatively long chords through the pith were assumed, measured and compared. Finally, the longest chord passing through the pith on the tree disc was obtained.

Perimeter measurement
In order to measure the length of every part of the outer contour of the tree disc, split by the longest chord through the pith, polylines were respectively made on the two parts, based on the principle of cyclotomy.    Sample images were shown in Figure 2. After sample images were imported into AutoCAD, the pith positions were marked manually. Then, as shown in Figure 3, many relatively long chords through the pith were assumed, measured and compared. Finally, the longest chord passing through the pith on the tree disc was obtained.

Perimeter measurement
In order to measure the length of every part of the outer contour of the tree disc, split by the longest chord through the pith, polylines were respectively made on the two parts, based on the principle of cyclotomy. 2 Perimeter measurement In order to measure the length of every part of the outer contour of the tree disc, split by the longest chord through the pith, polylines were respectively made on the two parts, based on the principle of cyclotomy. As shown in Figure 4, the outer-contour circumference, C was segmented into C 1 and C 2 . ∆C is defined as the length difference between C 1 and C 2 . ∆C of samples of different species are shown in Table 1.
Forests 2021, 12, x FOR PEER REVIEW 4 of 14 As shown in Figure 4, the outer-contour circumference, C was segmented into C1 and C2. ΔC is defined as the length difference between C1 and C2. ΔC of samples of different species are shown in Table 1.

t-test
Single sample t-test was used to analyze whether there was length difference between the two parts of the outer-contour circumference, split by the longest chord passing through the pith. The result of t-test shows that the significance of all species is less than 0.05, that is, the longest chord passing through the pith cannot bisect the outer-contour circumference without any error. In order to obtain the permissible error range, the range of ΔC and 95% confidence interval were analyzed, as shown in Table 2. In order to make an overall estimate of the error range, the upper limit of 95% confidence interval was employed during the error analysis, and the error E is defined as the formula (2). Error E of 12 samples in each species is shown in Table 3, and the 95% confidence interval of error E in each species is shown in Table 4.

t-test
Single sample t-test was used to analyze whether there was length difference between the two parts of the outer-contour circumference, split by the longest chord passing through the pith. The result of t-test shows that the significance of all species is less than 0.05, that is, the longest chord passing through the pith cannot bisect the outer-contour circumference without any error. In order to obtain the permissible error range, the range of ∆C and 95% confidence interval were analyzed, as shown in Table 2. In order to make an overall estimate of the error range, the upper limit of 95% confidence interval was employed during the error analysis, and the error E is defined as the Formula (2). Error E of 12 samples in each species is shown in Table 3, and the 95% confidence interval of error E in each species is shown in Table 4.  From Tables 3 and 4, it can be concluded that all the errors E of four different species are within 7%, lower than the permissible error range of 10% in this paper. Furthermore, single sample t-test with the test value of 5% was made to error E of each species samples. The result of t-test shows that the significance of all species is greater than 0.05, which means that there is no distinguished difference. Moreover, error E of samples of each species fluctuates around 5%, all of ∆C are within the permissible range of test error. Therefore, to an extent, the hypothesis about the geometric properties can be verified that the longest chord passing through the pith can bisect the outer-contour circumference of the tree disc within a permissible range of test error. In further research, more samples will be added to verify the geometric properties.

Data Acquisition
The data acquisition equipment is a 2D terrestrial laser scanner, SICK LMS141-15100 [28][29][30][31][32][33]. The laser lidar parameters are shown in Table 5. The experimental scene in the laboratory is shown in Figure 5. When using for standing trees in the field, the principle and the process are the same with that for tree discs. The tree disc was scanned at four points, including A, B, C, and D, with the same distance to O, the center of Square ABCD, as shown in Figure 6.   As seen in Figure 6, the position of the laser transmitter of LiDAR was regard the origin of the above-mentioned four coordinate systems, the direction orienting target disc as the direction of the positive axis X, and the direction of the positive was determined according to the right-hand rule. When the sample was scanned at A, B, C, and D, , , and , are defined as coordin the obtained data in the corresponding coordinate systems XAY, , . After four-time scanning of the tree disc, the point cloud data obtained by the L are shown in Figure 7. The blue dots represent the scanning data of point A, the re   As seen in Figure 6, the position of the laser transmitter of LiDAR was regarded as the origin of the above-mentioned four coordinate systems, the direction orienting to the target disc as the direction of the positive axis X, and the direction of the positive axis Y was determined according to the right-hand rule. When the sample was scanned at points A, B, C, and D, , , and , are defined as coordinates of the obtained data in the corresponding coordinate systems XAY, , , and . After four-time scanning of the tree disc, the point cloud data obtained by the LiDAR are shown in Figure 7. The blue dots represent the scanning data of point A, the red dots As seen in Figure 6, the position of the laser transmitter of LiDAR was regarded as the origin of the above-mentioned four coordinate systems, the direction orienting to the target disc as the direction of the positive axis X, and the direction of the positive axis Y was determined according to the right-hand rule. When the sample was scanned at points A, B, C, and D, (x Ai , y Ai ),(x Bi , y Bi ), (x Ci , y Ci ) and (x Di , y Di ) are defined as coordinates of the obtained data in the corresponding coordinate systems XAY, X BY , X CY and X DY .
After four-time scanning of the tree disc, the point cloud data obtained by the LiDAR are shown in Figure 7. The blue dots represent the scanning data of point A, the red dots represent the scanning data of point B, the green dots represent the scanning data of point C, and the yellow dots represent the scanning data of point D.
represent the scanning data of point B, the green dots represent the scanning data of point C, and the yellow dots represent the scanning data of point D.

Data Processing
The flow chart of data processing is listed in Figure 8.

Point Cloud Coordinate Transformation
From the spatial position relationship among the four Cartesian coordinate systems, it can be inferred that the transformation relationship between the different coordinate systems [34]. As shown in formulas (3)-(5), the transformation matrix, , , respectively means the matrix to transform the point cloud in the coordinate system , , to the coordinate system XAY.

Data Processing
The flow chart of data processing is listed in Figure 8.

Data Processing
The flow chart of data processing is listed in Figure 8.

Point Cloud Coordinate Transformation
From the spatial position relationship among the four Cartesian coordinate systems, it can be inferred that the transformation relationship between the different coordinate systems [34]. As shown in formulas (3)-(5), the transformation matrix, , , respectively means the matrix to transform the point cloud in the coordinate system , , to the coordinate system XAY.

Point Cloud Coordinate Transformation
From the spatial position relationship among the four Cartesian coordinate systems, it can be inferred that the transformation relationship between the different coordinate systems [34]. As shown in Formulas (3)-(5), the transformation matrix, T A B , T A C , T A D respectively means the matrix to transform the point cloud in the coordinate system X BY , X CY , X DY to the coordinate system XAY.
After transforming all the scanning data to the coordinate system XAY, the visualization results of four-time scanning data are displayed in Figure 9.
After transforming all the scanning data to the coordinate system XAY, the visualization results of four-time scanning data are displayed in Figure 9.

Point Cloud Denoising
Because there are lots of noise data in the raw data, denoising of point cloud in this paper was divided into two steps. First, the gross error was removed. According to the scanning principle of LiDAR, the measuring angle of the marginal part of LiDAR is bigger than the central one, so that there was a relatively large error in the marginal data compared with the central data, which led to the distortion of the marginal data. This paper filtered out the most marginal point cloud, and the first-step denoising result is shown in Figure 10. The second-step denoising made use of the distance threshold method [35][36][37]. That is to say, the data obtained from four-time scanning were grouped to denoise, according to the different scanning positions. Since the curve of each group of data is approximately arc-like, when the points on the arc satisfied formula (6) were retained. Otherwise, they were filtered out, and the second-step denoising result is presented in Figure 11.

Point Cloud Denoising
Because there are lots of noise data in the raw data, denoising of point cloud in this paper was divided into two steps. First, the gross error was removed. According to the scanning principle of LiDAR, the measuring angle of the marginal part of LiDAR is bigger than the central one, so that there was a relatively large error in the marginal data compared with the central data, which led to the distortion of the marginal data. This paper filtered out the most marginal point cloud, and the first-step denoising result is shown in Figure 10.
After transforming all the scanning data to the coordinate system XAY, the visualization results of four-time scanning data are displayed in Figure 9.

Point Cloud Denoising
Because there are lots of noise data in the raw data, denoising of point cloud in this paper was divided into two steps. First, the gross error was removed. According to the scanning principle of LiDAR, the measuring angle of the marginal part of LiDAR is bigger than the central one, so that there was a relatively large error in the marginal data compared with the central data, which led to the distortion of the marginal data. This paper filtered out the most marginal point cloud, and the first-step denoising result is shown in Figure 10. The second-step denoising made use of the distance threshold method [35][36][37]. That is to say, the data obtained from four-time scanning were grouped to denoise, according to the different scanning positions. Since the curve of each group of data is approximately arc-like, when the points on the arc satisfied formula (6) were retained. Otherwise, they were filtered out, and the second-step denoising result is presented in Figure 11. The second-step denoising made use of the distance threshold method [35][36][37]. That is to say, the data obtained from four-time scanning were grouped to denoise, according to the different scanning positions. Since the curve of each group of data is approximately arc-like, when the points on the arc satisfied Formula (6) were retained. Otherwise, they were filtered out, and the second-step denoising result is presented in Figure 11. where d i is defined as the distance between the data point and the center point O, d means the average distance of the group, and n is the number of collected point cloud.
where di is defined as the distance between the data point and the center point O, means the average distance of the group, and n is the number of collected point cloud. Figure 11. The result of the second-step denoising.

Point Cloud Sorting
The order of point cloud plays an important role in curve reconstruction [38]. The ant colony algorithm was employed to sort the denoised data. As seen in Figure 12a, the connected polygon of the unsorted point cloud is relatively unsmooth, and the connected sequence of several data points is obviously wrong. After sorting by the ant colony algorithm [39][40][41], it can be seen that the connected order of the point cloud is significantly improved in Figure 12b.

Point Cloud Fitting of Cross-Section Contour
In this paper, cubic spline curves were adopted to interpolate point cloud data of cross-section outline under the cyclic boundary condition [42][43][44]. The data visualization results after interpolation are shown in Figure 13. The blue-marked points are the original data, and the red-marked points are 200 interpolation points.
The point cloud data were connected to generate the outer-contour fitting curve of the cross-section, as shown in the blue line of Figure 13. Compared with the real contour

Point Cloud Sorting
The order of point cloud plays an important role in curve reconstruction [38]. The ant colony algorithm was employed to sort the denoised data. As seen in Figure 12a, the connected polygon of the unsorted point cloud is relatively unsmooth, and the connected sequence of several data points is obviously wrong. After sorting by the ant colony algorithm [39][40][41], it can be seen that the connected order of the point cloud is significantly improved in Figure 12b.
where di is defined as the distance between the data point and the center point O, means the average distance of the group, and n is the number of collected point cloud. Figure 11. The result of the second-step denoising.

Point Cloud Sorting
The order of point cloud plays an important role in curve reconstruction [38]. The ant colony algorithm was employed to sort the denoised data. As seen in Figure 12a, the connected polygon of the unsorted point cloud is relatively unsmooth, and the connected sequence of several data points is obviously wrong. After sorting by the ant colony algorithm [39][40][41], it can be seen that the connected order of the point cloud is significantly improved in Figure 12b.

Point Cloud Fitting of Cross-Section Contour
In this paper, cubic spline curves were adopted to interpolate point cloud data of cross-section outline under the cyclic boundary condition [42][43][44]. The data visualization results after interpolation are shown in Figure 13. The blue-marked points are the original data, and the red-marked points are 200 interpolation points.
The point cloud data were connected to generate the outer-contour fitting curve of the cross-section, as shown in the blue line of Figure 13. Compared with the real contour

Point Cloud Fitting of Cross-Section Contour
In this paper, cubic spline curves were adopted to interpolate point cloud data of cross-section outline under the cyclic boundary condition [42][43][44]. The data visualization results after interpolation are shown in Figure 13. The blue-marked points are the original data, and the red-marked points are 200 interpolation points.
of the disc, it can be found that both are basically coincident, indicating that the performance of interpolation and curve-fitting are good.

Finding the Longest Chord
As shown in formula (7), the chord lengths between two points were calculated and then the longest chord on the cross-section was obtained.
where , and , are any two points, and CLi is defined as the chord length between the corresponding two points. As seen in Figure 14, CLmax, with the two endpoints Nmax1 and Nmax2, is the longest chord, assumed to pass through the pith.

The Radial Location of the Pith
Supposed the fitting perimeter of the disc's cross-section is S, the perimeter of the two parts, segmented by the assumed longest chord, is respectively S1 and S2. The fitting curve is divided into several sub-intervals, composed of two adjacent points on the curve. As seen in Figure 15, the length of the curve in each sub-interval is ∆ and the length of the straight line is ∆ , and the differences of x-coordinate and y-coordinate at the endpoints are separately dx and dy. The length of the curve in each sub-interval is simplified as the length of the straight line between the two endpoints. It can be expressed as formulas (8) and (9). The point cloud data were connected to generate the outer-contour fitting curve of the cross-section, as shown in the blue line of Figure 13. Compared with the real contour of the disc, it can be found that both are basically coincident, indicating that the performance of interpolation and curve-fitting are good.

Finding the Longest Chord
As shown in Formula (7), the chord lengths between two points were calculated and then the longest chord on the cross-section was obtained.
where (x i , y i ) and x j , y j are any two points, and CL i is defined as the chord length between the corresponding two points. As seen in Figure 14, CL max , with the two endpoints N max1 and N max2 , is the longest chord, assumed to pass through the pith.

Finding the Longest Chord
As shown in formula (7), the chord lengths between two points were calculated and then the longest chord on the cross-section was obtained.
where , and , are any two points, and CLi is defined as the chord length between the corresponding two points. As seen in Figure 14, CLmax, with the two endpoints Nmax1 and Nmax2, is the longest chord, assumed to pass through the pith.

The Radial Location of the Pith
Supposed the fitting perimeter of the disc's cross-section is S, the perimeter of the two parts, segmented by the assumed longest chord, is respectively S1 and S2. The fitting curve is divided into several sub-intervals, composed of two adjacent points on the curve. As seen in Figure 15, the length of the curve in each sub-interval is ∆ and the length of the straight line is ∆ , and the differences of x-coordinate and y-coordinate at the endpoints are separately dx and dy. The length of the curve in each sub-interval is simplified as the length of the straight line between the two endpoints. It can be expressed as formulas (8) and (9).

The Radial Location of the Pith
Supposed the fitting perimeter of the disc's cross-section is S, the perimeter of the two parts, segmented by the assumed longest chord, is respectively S 1 and S 2 . The fitting curve is divided into several sub-intervals, composed of two adjacent points on the curve. As seen in Figure 15, the length of the curve in each sub-interval is ∆S i and the length of the straight line is ∆L i , and the differences of x-coordinate and y-coordinate at the endpoints are separately dx and dy. The length of the curve in each sub-interval is simplified as the length of the straight line between the two endpoints. It can be expressed as Formulas (8) and (9).
Forests 2021, 12, x FOR PEER REVIEW 11 of 14 ∆S i ≈∆L i = dx 2 + dy 2 (8) Based on the calculation results of S1 and S2, if the longest chord could not bisect the circumference of the cross-section, it was concluded that this chord does not pass through the pith. The second longest chord on the cross-section need to be searched until a longer chord conforming to the geometric properties, that is, bisecting the circumference, was found. Consequently, the radial direction of the pith was located.

Experiment and Analysis
In this paper, discs were employed to make location experiments out of cross-section because the visible pith made it easier to evaluate the location performance. Due to water loss over time, half of 10 tree discs had cracks as shown in Figure 16. In order to make sure that the above-mentioned geometric properties of tree pith apply to the cracked tree discs, error E of ΔC, the length difference between the two outer-contour circumferences, split by the longest chord passing through the pith, was calculated for the digital images of 10 tree discs. The results are shown in Table 6. In Table 6, all the errors of ΔC of 10 tree discs are lower than 5%, that is, the geometric properties of tree pith work for 10 tree discs. Therefore, 10 tree disc samples were tested outside of the cross-section for radial pith location based on LiDAR. The results are shown in Table 7. The error ECL of the two-part outer contour circumference, S1 and S2, is defined as formula (10).
According to the regulation for trees with DBH (diameter at breast height) greater than 20 cm in 'the Technical Regulations on Continuous Inventory of National Forest Resources' (State Forestry Administration, 2014), the DBH measurement error is required to be less than 1.5%. Therefore, the permissible error range in this paper is defined to be 1.5%.  Based on the calculation results of S 1 and S 2 , if the longest chord could not bisect the circumference of the cross-section, it was concluded that this chord does not pass through the pith. The second longest chord on the cross-section need to be searched until a longer chord conforming to the geometric properties, that is, bisecting the circumference, was found. Consequently, the radial direction of the pith was located.

Experiment and Analysis
In this paper, discs were employed to make location experiments out of cross-section because the visible pith made it easier to evaluate the location performance. Due to water loss over time, half of 10 tree discs had cracks as shown in Figure 16. In order to make sure that the above-mentioned geometric properties of tree pith apply to the cracked tree discs, error E of ∆C, the length difference between the two outer-contour circumferences, split by the longest chord passing through the pith, was calculated for the digital images of 10 tree discs. The results are shown in Table 6. CLmax of Sample 8 and 10 through the pith is respectively the second and third longest one of all the chords on the cross-section, and CLmax of all the other samples through the pith are the longest chord. As shown in Table 7, except for Sample 4 and 9, ECL errors of other samples are less than 1.5%, and the minimum error is only 0.17%. That is, for the other eight samples except Samples 4 and 9, the radial location of the pith can be achieved by finding the longest chord that bisects the circumference of the cross-section contour of the disc. On the cross-sections of Samples 4 and 9, the first five longest chords bisect the circumference of the cross-section outer contour of the disc with an error of more than 1.5%. A comparison was made between the experimental result and the ground truth of the sample discs, as shown in Figure 16. The red line means the chord passing through the pith obtained by the experiment, and the black dot was marked manually, representing the real pith of the disc. In Figure 16, it was concluded that all the CLmax obtained by this method passed through the pith O on the cross-section of the sample discs. In Sample 4 and 9, CLmax also passed through the real pith, although the bisectional effect of the circumference of the cross-section outer contour was not ideal.
By comparing the cross-section shapes of all the sample discs, it was found that the cross-section shapes of Sample 4 and 9 are more similar to circles, while other samples are  In Table 6, all the errors of ∆C of 10 tree discs are lower than 5%, that is, the geometric properties of tree pith work for 10 tree discs. Therefore, 10 tree disc samples were tested outside of the cross-section for radial pith location based on LiDAR. The results are shown in Table 7. The error E CL of the two-part outer contour circumference, S 1 and S 2 , is defined as Formula (10). According to the regulation for trees with DBH (diameter at breast height) greater than 20 cm in 'the Technical Regulations on Continuous Inventory of National Forest Resources' (State Forestry Administration, 2014), the DBH measurement error is required to be less than 1.5%. Therefore, the permissible error range in this paper is defined to be 1.5%.
CL max of Sample 8 and 10 through the pith is respectively the second and third longest one of all the chords on the cross-section, and CL max of all the other samples through the pith are the longest chord. As shown in Table 7, except for Sample 4 and 9, E CL errors of other samples are less than 1.5%, and the minimum error is only 0.17%. That is, for the other eight samples except Samples 4 and 9, the radial location of the pith can be achieved by finding the longest chord that bisects the circumference of the cross-section contour of the disc. On the cross-sections of Samples 4 and 9, the first five longest chords bisect the circumference of the cross-section outer contour of the disc with an error of more than 1.5%.
A comparison was made between the experimental result and the ground truth of the sample discs, as shown in Figure 16. The red line means the chord passing through the pith obtained by the experiment, and the black dot was marked manually, representing the real pith of the disc. In Figure 16, it was concluded that all the CL max obtained by this method passed through the pith O on the cross-section of the sample discs. In Sample 4 and 9, CL max also passed through the real pith, although the bisectional effect of the circumference of the cross-section outer contour was not ideal. By comparing the cross-section shapes of all the sample discs, it was found that the cross-section shapes of Sample 4 and 9 are more similar to circles, while other samples are more similar to ellipses. For the disc with the oval cross-section, the method proposed in this paper can be used to locate the radial direction of the pith in trees discs. However, for the disc with the circle-like cross-section shape, the longest chord on the cross-section would pass through the pith to implement radial location, and it is not necessary to verify whether this chord has the geometric property of bisecting the cross-section perimeter.

Conclusions
In this paper, the geometric property of the pith was first proposed and verified, that is, the longest chord passing through the pith, can bisect the circumference of cross-section outline of the disc within a certain error range. On this basis, a radial location method of the pith based on terrestrial laser scanning was put forward. Two-dimensional LiDAR was used to scan tree discs to obtain data. Then, coordinate transformation, denoising, sorting, and interpolation of the point cloud were carried out, and the longest chords bisecting the cross-section perimeter were obtained to locate the radial direction of the pith. Through experimental verification and comparative analysis on tree discs, for the tree disc whose cross-section shape is close to ellipse, the error of the method proposed in this paper is less than 1.5%.