# Hardwood Species Classification with Hyperspectral Microscopic Images

Authors: Peng Zhao; Cheng-Kun Wang
Journal: Journal of Spectroscopy
 (2019)

---

## Abstract

We propose a hardwood species identification method based on wood hyperspectral microscopic images. A SOC710VP hyperspectral stereomicroscope was used to acquire microscopic images of a hardwood cross section. In these microscopic images, each part’s spectral features are discussed. We found that the spectral divisibility of wood vessels’ peripheral and central regions in the hyperspectral microscopic images can be used for hardwood species recognition. Mathematical morphological operation and K-L divergence were used to extract spectral features at the wood vessels’ peripheral regions and central regions, respectively. By comparing wood vessels’ spectral similarity across wood species samples, we found that wood vessels’ peripheral spectral divisibility is larger than its central. Finally, the spectral information from randomly selected regions of interest (i.e., ROI) and that of wood vessels’ peripheral and central regions have value as a classification basis. In our hardwood species classification experiments, three dimensionality reduction algorithms, principal component analysis (PCA), kernel principal component analysis (KPCA), and multidimensional scaling (MDS), and the three classifiers, BP neural network, support vector machine (SVM), and Mahalanobis distance (MD), are combined to perform hardwood species classification work. Experimental results indicate that the best recognition effect can be achieved at the peripheral region of wood vessels using PCA or MDS with the MD algorithm.


---

## Body

## 1. Introduction

Wood is an important renewable natural resource. Different tree speciesresult in many different wood products circulating in the market. In particular, many hardwood species are widely used to make furniture, and they cannot be identified by visual observation or odor discrimination. Hardwood misclassification may result in significant economic losses due to the differences in quality and price of different hardwood species.With the rapid development of machine vision, it has become possible to identify the processed hardwood species by using intelligent methods. Pan and Kudo [1] and Yusof et al. [2] used microscopic images of wood cross sections to discuss the characteristics of different hardwood species. They showed that different hardwood species have differences in their wood’s vessel structure. Furthermore, some textural or statistical properties for wood vessels were extracted from the 2D wood images to perform hardwood species classification [3, 4]. Macroscopic wood features are also used in wood species recognition. For example, Hu et al. used a deep learning scheme for lumber grade classification in terms of lumber defects, textures, and species [5].In recent years, spectral analysis schemes have been proposed for wood species classification [6–8]. These schemes are in essence spectrum-based classifications which mainly use 1D spectral curves. They usually have low computational complexity with fast speed due to the 1D data but may suffer from spectral noise and disturbance from the surrounding site. Therefore, wood species recognition with 3D hyperspectral imaging which includes both visual and spectral information of wood tissues may be a good choice. In fact, 3D hyperspectral imaging schemes have been used in agricultural engineering for food quality evaluation [9] and crop content predictions [10, 11].In this paper, the microscopic hyperspectral images of hardwood cross sections were recorded using a SOC710VP hyperspectral stereomicroscope to investigate the possibility of hardwood species identifications. The spectral features of different wood tissues such as wood vessel, wood ray, and background are discussed. It was found that spectral information from the peripheral and central regions of wood vessels can be used for hardwood species identification.


## 2. Materials and Methods

### 2.1. Wood Samples

Six hardwood species were studied (Table1). One thousand four hundred and forty wood samples were prepared from six species, and each sample was cut into a cube with length, width, and height of 2 cm, 2 cm, and 3 cm, respectively. Nine hundred samples were randomly selected as the training set, and the other 540 samples are used as the test set. Each sample’s cross sections (2cm×2cm) were cut using a disk saw to make the two surfaces as smooth as possible, while the radial sections and tangential sections (3cm×2cm) were cut using an electric plane. The wood sample’s cross sections were used for hyperspectral imaging. In these cross sections, the wood growth rings were as close to parallel with the section’s two edges as possible. Quercus rubra and Fraxinus americana had similar colors, while cross sections for the other four species had different colors.Table 1
Species studied.
Number
English name
Scientific name
1
Red oak
Quercus rubra
2
Ash
Fraxinus americana
3
Merbau
Intsia bijuga
4
Golden silk teak
Tectona grandis
5
Rose wood
Guibourtia tessmannii
6
Limba
Terminalia neotaliala capuronThe intraclass differences of wood samples are caused by the tree’s age, cutting position, and place of growth. For each wood species, our wood blocks are taken from different individual trees at random position in the cross section of the timber to include these intraclass differences. In this way, the reliability of our proposed scheme can be ensured.


### 2.2. Spectral Data Acquisition

The spectral acquisition instrument used for this research was a SOC710VP hyperspectral stereomicroscope, with 128 spectral bands and a spectral range of 372–1038 nm. To improve spectral acquisition, all wood sample surfaces were flat, without burrs, and the spectral acquisition environment was consistent and stable. The dimension of captured microscopic hyperspectral images was520×696×128. To eliminate the influence of dark current, it was also necessary to perform black and white correction processing on the detected samples. The wood image’s enlargement factor was 45.To control costs, we chose the visible/near infrared hyperspectral stereomicroscope as the experimental instrument. Spectral data in the visible region may vary because wood sample color may change with environmental variation. Therefore, we maintained uniform environmental conditions during the spectral experiments and wood sample preservation. In addition, 780–900 nm is the near infrared band, which best characterizes the different wood properties. A stereomicroscope was used to capture wood microscopic images conveniently for each wood sample. The wood slice is not used as its production is complex and time consuming.


### 2.3. Spectral Analysis of Different Tissues of Hardwood Microscopic Hyperspectral Images

The microscopic hyperspectral images of cross sections mainly include wood vessels, wood rays, and axial parenchyma, as shown in Figure1. The closed holes are wood vessels which are important features for species recognition [12, 13]. The dense, thin, approximately parallel lines are called wood rays, while the thick lines approximately perpendicular to the wood rays are the axial parenchyma.Figure 1
Microscopic hyperspectral images of cross sections of different hardwood species. (a)Quercus rubra. (b) Fraxinus americana. (c) Intsia bijuga.
(a)
(b)
(c)Hyperspectral images are 3D cubes which have both spatial information and spectral information, and the spectral curves of different wood tissues for the same hardwood species are different. For example, Figure2 illustrates the spectral curves of Intsia bijuga in the central and peripheral regions of vessels, in the regions of wood rays, axial parenchyma, and the background. Since vessels are the most important features in hardwood species identification [1, 3, 4, 12], the spectral information in the vessels’ central and peripheral regions is extracted and discussed.Figure 2
Spectral curves for different tissues ofIntsia bijuga: (a) the 5 cross marks indicate different tissues; (b) corresponding spectral curves for each tissue.
(a)
(b)


### 2.4. Spectral Analysis in Vessel Peripheral Regions with Mathematical Morphology

Mathematical morphology is widely used in image processing, which can be used to extract useful image components to enhance ROI and suppress background regions. We denote the original image, background image, and ROI image asHx,y, Bx,y, and Fx,y, respectively. The structural element with a specified shape is denoted as A. The background image and ROI image can be expressed as follows:(1)Bx,y=Hx,y○Ax,y,Fx,y=Hx,y−Bx,y,where ○ represents the morphological opening operation.In order to eliminate noise, a corrosion operation is required for the obtainedFx,y. Then, the binary image is formed by using threshold processing to Fx,y. The gray-level image for each band of hyperspectral image is processed using the abovementioned procedures to obtain the 128 binary images which form a cube of M×N×128 dimensions (i.e., MandN are the width and height of hyperspectral image).This cube can be regarded asM×N vectors of 128 dimensions, and the element in every vector is either 0 or 1. The number of element 1 in every vector is registered so as to form a matrix D of M×N dimensions. A threshold k is defined to convert the matrix D into a binary image E as follows, which is illustrated in Figure 3.(2)Ex,y=1,ifDx,y≥k,0,ifDx,y<k.Figure 3
The wood vessels’ peripheral region extraction with mathematical morphology: (a) the hyperspectral 56th band image; (b) the binary image cube of M×N×128 dimensions; (c) the binary image E.
(a)
(b)
(c)For the binary imageE, the spectra are extracted whose corresponding pixel value is 1 so as to form a new matrix X of C×128 dimensions, while C represents the number of element 1 in binary image E.The abovementioned image and spectral processing operations were performed to a ROI in a hyperspectral image of one sample. In this work, the size of ROI is set as100×100. In a wood hyperspectral image of 520×696×128 dimensions, the initial ROI is placed in the upper-left position in the image to obtain the computed matrix X. Then, the ROI moves gradually from left to right and from top to bottom to traverse the whole hyperspectral image to obtain multiple matrices Xi,i=1,2,…,n. These matrices Xi,i=1,2…n are connected by columns to obtain a matrix Y with a dimension of C1+⋯+Cn×128, and the mean value for each column is computed to form a new vector y of 128 dimensions. For a hardwood species with m hyperspectral image specimens, we may obtain multiple vectors yj,j=1,2,…,m and calculate the mean vector y˜. As an example, Figure 4 illustrates the mean spectral vector y˜ for the six hardwood species in our measurements in the wood vessels’ peripheral regions. The training samples are used to calculate y˜ with the ROI of 100×100 and the structural element A of disk shape with the size r. In Figure 4, we can see that the mean spectral curves for training sets in wood vessels’ peripheral regions are different for the six hardwood species. Therefore, this spectral difference can be used to identify hardwood species. Since different hardwood species have different vessel sizes, the structural element A with different size r is used for the 6 hardwood species. This is explained in detail in Section 3.4.Figure 4
Spectral difference in wood vessels’ peripheral region for six hardwood species.


### 2.5. Spectral Analysis in Wood Vessels’ Central Regions with K-L Divergence

K-L divergence is widely used in information theory as a signal similarity measuring tool [14]. We suppose the probability distributions of two discrete random signals as X=x1x2…xnT and Y=y1y2…ynT, respectively. The K-L divergence for Y to X is then defined as follows:(3)DKLXY=∑i=1nxilogaxiyi,∑i=1nxi=∑i=1nyi=1.Therefore, the difference between two band images of a hyperspectral image cube can be measured using the K-L divergence equation (3).Suppose a hyperspectral image cubeX=x1,x2,…,xLT∈RL×N, where xi=xi1,xi2,…,xiN and L, N represent the number of spectral bands and pixels of a 2D band image, respectively. After normalization of data in X, the K-L divergence of the jth band image relative to the ith band image can be defined as follows:(4)DKLxixj=∑k=1Nxiklogxikxjk.In this way, we can calculate the K-L divergence of two band images to obtain a symmetric matrixD. The ith row of this matrix is defined as ai=Di,1…Di,L, and the mean value of its L elements represents, to some extent, the information entropy of the ith band image. Therefore, we may use the 10 spectral bands with the 10 largest entropies to form an image cube X˜=xλ1,xλ2,…,xλ10T∈R10×N (Figure 5(a)) and its complementary image cube X⌣ (Figure 5(b)). This cube X⌣ is classified into 4 parts by use of a K-means clustering algorithm so that the spatial pixels in the original wood hyperspectral image cube are classified into 4 parts, as illustrated in Figures 5(c)–5(f).Figure 5
Spatial clustering results with spectral information andK-means algorithm: (a) image cube with selected 10 bands; (b) the complement image cube; (c–f) the classified 4 parts.
(a)
(b)
(c)
(d)
(e)
(f)As shown in Figure5, the white part in Figure 5(d) is the wood vessels’ central region, which is segmented accurately. Figures 5(c), 5(e) and 5(f) indicate the wood vessels’ peripheral and background regions. In fact, the wood vessels’ peripheral region is clustered separately into different parts, as illustrated in Figures 5(c), 5(e) and 5(f). Therefore, this classification is not suitable for accurately extracting the wood vessels’ peripheral region. In order to more accurately segment the vessels’ central region, it is necessary to perform a morphological corrosion operation on the white sections in Figure 5(d).The clustering process in Figure5 has two problems. Firstly, selecting the ROI randomly in a hyperspectral image does not ensure location of the wood vessels’ central region. Secondly, as can be seen from Figure 5, the K-means cluster results in 4 parts (i.e., Figures 5(c)–5(f)), but only Figure 5(d) represents the wood vessels’ central region.To solve these problems, we propose an automatic identification scheme for the wood vessels’ central region. For a ROIW of size S×S, the binary images produced by K-means cluster are B1,B2,B3,B4. The selected 10 spectral bands obtained from K-L divergence are λ1,λ2,…,λ10, and the corresponding binary images for these 10 band images are denoted as Xn,n=1,2,3,…,10. Here CpX,p=0,1 represents the occurrence number of p in matrix X, and CX,W represents the number of identical elements in the same position of matrices X and W. The 4 measurement parameters are defined as follows:(5)b1k=1/10∑i=1i=10C0XiC0Bk,k∈1,2,3,4,(6)b2k=1/10∑i=1i=10C1XiC1Bk,k∈1,2,3,4,(7)b3k=1/10∑i=1i=10CXi,BkS×S,k∈1,2,3,4,(8)bk=13b1k+b2k+b3k.Herex is defined as follows:(9)x=x,x≤1,1x,x>1.We calculatedbk for the binary image Bk,k∈1,2,3,4 with equation (8) and define g=maxbk and a threshold T. When g≥T, the corresponding binary image Bk includes the wood vessels’ central region. Using this approach, the ROI moves gradually from left to right and from top to bottom with a fixed step length S to traverse the whole hyperspectral image. In this way, we locate all binary images Ki,i=1,2,…,l and include the wood vessels’ central region. These image matrices Ki,i=1,2,…,l are connected by columns to obtain a matrix Z with a dimension of C1+⋯+Cl×128, and the mean value of each column is computed to form a new vector z of 128 dimensions. For a hardwood species with m hyperspectral image specimens, we derive multiple vectors zj,j=1,2,…,m and calculate the mean vector z˜. Figure 6 illustrates the mean spectral vector z˜ for the six hardwood species in our studies of the wood vessels’ central regions. The training samples are used to calculate z˜ with a ROI size of 100×100.Figure 6
Spectral difference in wood vessels’ central region for six hardwood species.


## 2.1. Wood Samples

Six hardwood species were studied (Table1). One thousand four hundred and forty wood samples were prepared from six species, and each sample was cut into a cube with length, width, and height of 2 cm, 2 cm, and 3 cm, respectively. Nine hundred samples were randomly selected as the training set, and the other 540 samples are used as the test set. Each sample’s cross sections (2cm×2cm) were cut using a disk saw to make the two surfaces as smooth as possible, while the radial sections and tangential sections (3cm×2cm) were cut using an electric plane. The wood sample’s cross sections were used for hyperspectral imaging. In these cross sections, the wood growth rings were as close to parallel with the section’s two edges as possible. Quercus rubra and Fraxinus americana had similar colors, while cross sections for the other four species had different colors.Table 1
Species studied.
Number
English name
Scientific name
1
Red oak
Quercus rubra
2
Ash
Fraxinus americana
3
Merbau
Intsia bijuga
4
Golden silk teak
Tectona grandis
5
Rose wood
Guibourtia tessmannii
6
Limba
Terminalia neotaliala capuronThe intraclass differences of wood samples are caused by the tree’s age, cutting position, and place of growth. For each wood species, our wood blocks are taken from different individual trees at random position in the cross section of the timber to include these intraclass differences. In this way, the reliability of our proposed scheme can be ensured.


## 2.2. Spectral Data Acquisition

The spectral acquisition instrument used for this research was a SOC710VP hyperspectral stereomicroscope, with 128 spectral bands and a spectral range of 372–1038 nm. To improve spectral acquisition, all wood sample surfaces were flat, without burrs, and the spectral acquisition environment was consistent and stable. The dimension of captured microscopic hyperspectral images was520×696×128. To eliminate the influence of dark current, it was also necessary to perform black and white correction processing on the detected samples. The wood image’s enlargement factor was 45.To control costs, we chose the visible/near infrared hyperspectral stereomicroscope as the experimental instrument. Spectral data in the visible region may vary because wood sample color may change with environmental variation. Therefore, we maintained uniform environmental conditions during the spectral experiments and wood sample preservation. In addition, 780–900 nm is the near infrared band, which best characterizes the different wood properties. A stereomicroscope was used to capture wood microscopic images conveniently for each wood sample. The wood slice is not used as its production is complex and time consuming.


## 2.3. Spectral Analysis of Different Tissues of Hardwood Microscopic Hyperspectral Images

The microscopic hyperspectral images of cross sections mainly include wood vessels, wood rays, and axial parenchyma, as shown in Figure1. The closed holes are wood vessels which are important features for species recognition [12, 13]. The dense, thin, approximately parallel lines are called wood rays, while the thick lines approximately perpendicular to the wood rays are the axial parenchyma.Figure 1
Microscopic hyperspectral images of cross sections of different hardwood species. (a)Quercus rubra. (b) Fraxinus americana. (c) Intsia bijuga.
(a)
(b)
(c)Hyperspectral images are 3D cubes which have both spatial information and spectral information, and the spectral curves of different wood tissues for the same hardwood species are different. For example, Figure2 illustrates the spectral curves of Intsia bijuga in the central and peripheral regions of vessels, in the regions of wood rays, axial parenchyma, and the background. Since vessels are the most important features in hardwood species identification [1, 3, 4, 12], the spectral information in the vessels’ central and peripheral regions is extracted and discussed.Figure 2
Spectral curves for different tissues ofIntsia bijuga: (a) the 5 cross marks indicate different tissues; (b) corresponding spectral curves for each tissue.
(a)
(b)


## 2.4. Spectral Analysis in Vessel Peripheral Regions with Mathematical Morphology

Mathematical morphology is widely used in image processing, which can be used to extract useful image components to enhance ROI and suppress background regions. We denote the original image, background image, and ROI image asHx,y, Bx,y, and Fx,y, respectively. The structural element with a specified shape is denoted as A. The background image and ROI image can be expressed as follows:(1)Bx,y=Hx,y○Ax,y,Fx,y=Hx,y−Bx,y,where ○ represents the morphological opening operation.In order to eliminate noise, a corrosion operation is required for the obtainedFx,y. Then, the binary image is formed by using threshold processing to Fx,y. The gray-level image for each band of hyperspectral image is processed using the abovementioned procedures to obtain the 128 binary images which form a cube of M×N×128 dimensions (i.e., MandN are the width and height of hyperspectral image).This cube can be regarded asM×N vectors of 128 dimensions, and the element in every vector is either 0 or 1. The number of element 1 in every vector is registered so as to form a matrix D of M×N dimensions. A threshold k is defined to convert the matrix D into a binary image E as follows, which is illustrated in Figure 3.(2)Ex,y=1,ifDx,y≥k,0,ifDx,y<k.Figure 3
The wood vessels’ peripheral region extraction with mathematical morphology: (a) the hyperspectral 56th band image; (b) the binary image cube of M×N×128 dimensions; (c) the binary image E.
(a)
(b)
(c)For the binary imageE, the spectra are extracted whose corresponding pixel value is 1 so as to form a new matrix X of C×128 dimensions, while C represents the number of element 1 in binary image E.The abovementioned image and spectral processing operations were performed to a ROI in a hyperspectral image of one sample. In this work, the size of ROI is set as100×100. In a wood hyperspectral image of 520×696×128 dimensions, the initial ROI is placed in the upper-left position in the image to obtain the computed matrix X. Then, the ROI moves gradually from left to right and from top to bottom to traverse the whole hyperspectral image to obtain multiple matrices Xi,i=1,2,…,n. These matrices Xi,i=1,2…n are connected by columns to obtain a matrix Y with a dimension of C1+⋯+Cn×128, and the mean value for each column is computed to form a new vector y of 128 dimensions. For a hardwood species with m hyperspectral image specimens, we may obtain multiple vectors yj,j=1,2,…,m and calculate the mean vector y˜. As an example, Figure 4 illustrates the mean spectral vector y˜ for the six hardwood species in our measurements in the wood vessels’ peripheral regions. The training samples are used to calculate y˜ with the ROI of 100×100 and the structural element A of disk shape with the size r. In Figure 4, we can see that the mean spectral curves for training sets in wood vessels’ peripheral regions are different for the six hardwood species. Therefore, this spectral difference can be used to identify hardwood species. Since different hardwood species have different vessel sizes, the structural element A with different size r is used for the 6 hardwood species. This is explained in detail in Section 3.4.Figure 4
Spectral difference in wood vessels’ peripheral region for six hardwood species.


## 2.5. Spectral Analysis in Wood Vessels’ Central Regions with K-L Divergence

K-L divergence is widely used in information theory as a signal similarity measuring tool [14]. We suppose the probability distributions of two discrete random signals as X=x1x2…xnT and Y=y1y2…ynT, respectively. The K-L divergence for Y to X is then defined as follows:(3)DKLXY=∑i=1nxilogaxiyi,∑i=1nxi=∑i=1nyi=1.Therefore, the difference between two band images of a hyperspectral image cube can be measured using the K-L divergence equation (3).Suppose a hyperspectral image cubeX=x1,x2,…,xLT∈RL×N, where xi=xi1,xi2,…,xiN and L, N represent the number of spectral bands and pixels of a 2D band image, respectively. After normalization of data in X, the K-L divergence of the jth band image relative to the ith band image can be defined as follows:(4)DKLxixj=∑k=1Nxiklogxikxjk.In this way, we can calculate the K-L divergence of two band images to obtain a symmetric matrixD. The ith row of this matrix is defined as ai=Di,1…Di,L, and the mean value of its L elements represents, to some extent, the information entropy of the ith band image. Therefore, we may use the 10 spectral bands with the 10 largest entropies to form an image cube X˜=xλ1,xλ2,…,xλ10T∈R10×N (Figure 5(a)) and its complementary image cube X⌣ (Figure 5(b)). This cube X⌣ is classified into 4 parts by use of a K-means clustering algorithm so that the spatial pixels in the original wood hyperspectral image cube are classified into 4 parts, as illustrated in Figures 5(c)–5(f).Figure 5
Spatial clustering results with spectral information andK-means algorithm: (a) image cube with selected 10 bands; (b) the complement image cube; (c–f) the classified 4 parts.
(a)
(b)
(c)
(d)
(e)
(f)As shown in Figure5, the white part in Figure 5(d) is the wood vessels’ central region, which is segmented accurately. Figures 5(c), 5(e) and 5(f) indicate the wood vessels’ peripheral and background regions. In fact, the wood vessels’ peripheral region is clustered separately into different parts, as illustrated in Figures 5(c), 5(e) and 5(f). Therefore, this classification is not suitable for accurately extracting the wood vessels’ peripheral region. In order to more accurately segment the vessels’ central region, it is necessary to perform a morphological corrosion operation on the white sections in Figure 5(d).The clustering process in Figure5 has two problems. Firstly, selecting the ROI randomly in a hyperspectral image does not ensure location of the wood vessels’ central region. Secondly, as can be seen from Figure 5, the K-means cluster results in 4 parts (i.e., Figures 5(c)–5(f)), but only Figure 5(d) represents the wood vessels’ central region.To solve these problems, we propose an automatic identification scheme for the wood vessels’ central region. For a ROIW of size S×S, the binary images produced by K-means cluster are B1,B2,B3,B4. The selected 10 spectral bands obtained from K-L divergence are λ1,λ2,…,λ10, and the corresponding binary images for these 10 band images are denoted as Xn,n=1,2,3,…,10. Here CpX,p=0,1 represents the occurrence number of p in matrix X, and CX,W represents the number of identical elements in the same position of matrices X and W. The 4 measurement parameters are defined as follows:(5)b1k=1/10∑i=1i=10C0XiC0Bk,k∈1,2,3,4,(6)b2k=1/10∑i=1i=10C1XiC1Bk,k∈1,2,3,4,(7)b3k=1/10∑i=1i=10CXi,BkS×S,k∈1,2,3,4,(8)bk=13b1k+b2k+b3k.Herex is defined as follows:(9)x=x,x≤1,1x,x>1.We calculatedbk for the binary image Bk,k∈1,2,3,4 with equation (8) and define g=maxbk and a threshold T. When g≥T, the corresponding binary image Bk includes the wood vessels’ central region. Using this approach, the ROI moves gradually from left to right and from top to bottom with a fixed step length S to traverse the whole hyperspectral image. In this way, we locate all binary images Ki,i=1,2,…,l and include the wood vessels’ central region. These image matrices Ki,i=1,2,…,l are connected by columns to obtain a matrix Z with a dimension of C1+⋯+Cl×128, and the mean value of each column is computed to form a new vector z of 128 dimensions. For a hardwood species with m hyperspectral image specimens, we derive multiple vectors zj,j=1,2,…,m and calculate the mean vector z˜. Figure 6 illustrates the mean spectral vector z˜ for the six hardwood species in our studies of the wood vessels’ central regions. The training samples are used to calculate z˜ with a ROI size of 100×100.Figure 6
Spectral difference in wood vessels’ central region for six hardwood species.


## 3. Results and Discussion

### 3.1. Spectral Difference Analysis for Wood Vessels’ Peripheral Region, Central Region, and Randomly Selected ROI

We may use spectral properties of a specific wood region or tissue to guide wood species identification using both interclass and intraclass differences. It is obvious that large interclass differences and small intraclass differences will produce more reliable results. These two differences can be represented by the spectral curve’s interclass and intraclass distances. Therefore, here, we use the divergence matrix to indicate the spectral curve’s interclass and intraclass differences. We define the intraclass divergence matrix for thejth wood species as follows:(10)Sj=1nj∑k=1njXkj−mjXkj−mjT.Herej=1,2,…,cc=6; Xkj represents the kth sample for the jth species, mj represents the mean value of all samples for the jth species, and nj represents the number of all samples for the jth species. The overall intraclass divergence matrix can be defined as follows:(11)Sw=∑j=1cPjSj,where Pj represents the prior probability for the jth wood species, and it can be estimated by Pj=nj/n, where n represents the overall number of all species samples. Similarly, the interclass divergence matrix is defined as follows:(12)Sb=∑j=1cPjmj−mmj−mT,where m represents the mean value of all species samples. The overall divergence matrix is then defined as follows:(13)St=Sb+Sw.Finally, the pattern divisibility rule is defined as follows:(14)J=trSbtrSw,where trX represents the trace of matrix X. It is obvious that a higher J value will produce a more reliable identification.In terms of a specific wood region or tissue, we considered the wood vessels’ peripheral region, central region, and randomly selected ROI. Spectral extractions in wood vessels’ peripheral and central regions have been discussed in Sections2.4 and 2.5. The size of the randomly selected ROI was 300×300; the mean vector of all spectral vectors in this ROI is the final spectral extraction for a hyperspectral cube sample.The smoothing procedure for the spectral preprocessing and the SNV calibration procedure are usually required before spectral extraction and pattern recognition as they cancel the noise disturbance and waveform offset. This preprocessing process is illustrated in Figure7.Figure 7
Spectral preprocessing procedure: (a) original spectral curve; (b) smoothed spectral curve; (c) spectral curve with SNV processing.
(a)
(b)
(c)Table2 gives the J values for the wood vessels’ peripheral region, central region, and randomly selected ROI in the training samples. It indicates that the spectral divisibility for the wood vessels’ peripheral region was most useful, while that for the wood vessels’ central region was least useful. Since the wood vessels’ central region is hollow, with some resin and probably some filling material, its spectral curve is variable and thus has a large intraclass difference with a small J value.Table 2
J values for different regions of the vessels.
Item
Vessels’ central region
Vessels’ peripheral region
Random ROI
J
1.6506
9.0970
6.6387Finally, the time required was compared for spectral extractions in the wood vessels’ peripheral region, central region, and random ROI. The computer configuration used was as follows: CPU Intel Core I7-6700, internal memory 8G, hard disk 1T, display card AMD Radeon R7200 Series. In these three spectral extractions, one hyperspectral image was cut into 30 ROI of100×100. Time required for processing one ROI is given in Table 3. We can see that the processing speed of random ROI was the fastest while that for the vessels’ central region was the slowest.Table 3
Time requirement comparisons for 3 spectral extractions.
Method
Process
Computing time (s)
Vessels’ peripheral region
G
x
,
y
0.6095
D
x
,
y
0.0056
E
x
,
y
0.0008
Total
0.6199
Vessels’ central region
D
0.7714
K-means cluster
0.2252
Central region selection
0.0087
Total
1.0181
Random ROI
Total
0.0874


### 3.2. Comparisons of Identification Accuracy

Before sending the extracted spectral vectors into classifiers, spectral dimensionality reduction was usually required. In this work, three algorithms, PCA, KPCA, and MDS, were used to reduce dimensionality. In spectral classification, three classifiers, BP neural networks, SVM, and MD, were used for identification comparisons.Firstly, we used the MD classifier to calculate the MD between a test sample and the training sample sets of the 6 hardwood species. The usual nearest neighbor classification was used to judge the test sample’s wood species. The classification results are illustrated in Figure8. We can see that the classification correct ratio (i.e., CCR) was influenced by spectral dimensionality reduction. For instance, if MDS was used, the CCR stabilized only when the spectral dimensionality was larger than 7. On the other hand, if PCA or KPCA was used, the CCR stabilized only when the spectral dimensionality was larger than 3. Moreover, we can see that the CCR was also influenced by the wood region of the spectral collection. In fact, the CCR of the wood vessels’ peripheral region was the highest, approaching 100%, while the CCR of the wood vessels’ central region was the lowest, and its mean value was approximately 70%–80%. The CCR of random ROI fell between the peripheral and central vessel regions and was approximately 90%.Figure 8
Spectral classification comparisons with MD classifier.The classification accuracy was the calculated CCR, which was the overall classification ratio (%) for the 540 test set samples. This classification accuracy is illustrated in Figures8–11 and Table 4.Figure 9
Identification accuracy of BP network for 25 training times (i.e., 25 networks).Figure 10
Spectral classification comparisons with BP neural network classifier.Figure 11
Spectral classification comparisons with SVM classifier.Table 4
Best classification results with different classifiers and spectral dimensionality reductions.
Dimensionality reductions
PCA
MDS
KPCA
Classifiers
Spectral region
CCR (%)
Spectral dimension
CCR (%)
Spectral dimension
CCR (%)
Spectral dimension
BP
Central vessel region
81.78
13
80.42
13
79.4
13
Random ROI
97.80
9
97.57
9
94.25
9
Peripheral vessel region
98.67
7
98.14
8
97.65
9
MD
Central vessel region
81.29
14
78.14
11
66.48
12
Random ROI
88.88
13
91.11
11
94.44
15
Peripheral vessel region
100
8
100
9
98.70
7
SVM
Central vessel region
83.70
15
83.33
14
83.33
13
Random ROI
97.78
7
97.78
8
94.62
12
Peripheral vessel region
99.07
8
99.07
7
98.88
7Secondly, we used the BP neural network classifier. Since the threshold and weight for the network initialization were selected randomly, the classification accuracy of every trained network may fluctuate slightly. Figure9 illustrates the test CCR for the 25 trained networks, while the spectral dimensionality was reduced to 5D by PCA. We observed that the overall CCR was stable, though there were slight differences between each CCR across the 25 networks. Therefore, the mean CCR for 25 networks was calculated and used as the final CCR. Figure 10 illustrates the mean CCR for different spectral dimensionality reductions and different wood regions of spectral collection. We can see that CCRs from vessel central regions were lowest, with mean values of approximately 70%–80%, while those from peripheral regions and random ROI were almost identical, approaching 100%.Thirdly, the SVM classifier was used. The predictive power of SVM depends on the kernel function used. The radial basis kernel function was usually used. When determining the penalty factorC and γ value, the grid search method can be used to find the best parameters. Figure 11 illustrates the identification comparisons with the best parameters C and γ. We can see again that the CCR of central regions was lowest, with a mean of approximately 80%, while those of peripheral regions and ROI were almost identical, approaching 100%. The CCR of random ROI decreased with the application of KPCA.The best identifications with different classifiers and spectral dimensionality reductions are given in Table4.Figures8, 10, and 11 and Table 4 indicated that the CCR of vessels’ central regions was the lowest in all cases and that the CCR of vessels’ peripheral regions was slightly higher than that of random ROI in the case of BP network and SVM. In the case of MD classifier, the CCR of vessels’ peripheral regions was much higher than that of random ROI. Moreover, we obtained good CCR of peripheral regions with a low spectral dimension. For example, the best CCR (100% of vessels’ peripheral region) was achieved using PCA and MD classifier with a spectral dimension of 8 or by using MDS and MD classifier with a spectral dimension of 9.Concerning the spectral dimensionality reductions, good results were obtained using PCA and KPCA with a low spectral dimension. However, if MDS was used, the CCR stabilized only when the spectral dimensionality was larger than 7. In general, PCA was the best choice in terms of identification accuracy and spectral dimension.


### 3.3. Comparisons of Identification Robustness

Wood hyperspectral image acquisition may be influenced by external environmental factors such as temperature, humidity, and illumination. In this section, a simulation experiment on noisy wood hyperspectral image classification was performed. Since the training set of wood images can be acquired in advance in a controlled environment, we only added white noises to the test set. The signal-to-noise ratio (SNR) was used to measure the test set’s image quality.The robustness of our three spectral extraction schemes (i.e. vessels’ central, vessels’ peripheral, random ROI) was compared using an error rate defined in terms of SNR.(15)Error rate=CCRideal−CCRnoisy,where CCRideal and CCRnoisy refer to the classification accuracy for ideal test images and noisy test images, respectively. The falling trend of “error rate” (Figure 12) indicates that the faster the curve falls, the more robust the corresponding spectral extraction scheme. Therefore, we conclude that in the case of BP and MD classifiers, the vessels’ peripheral extraction scheme was the most robust, while for SVM, the vessels’ peripheral scheme was slightly more robust than random ROI and much more robust than the vessels’ central scheme.Figure 12
Different SNR with error rates: (a) SVM; (b) BP; (c) MD.
(a)
(b)
(c)


### 3.4. Discussions of Parameters

In the wood vessels’ peripheral region extraction, the two key parameters are the sizer of structure element A and the threshold k in equation (2). Here we explain the effect of these two parameters and give the best values of r and k for each species.Figure13 illustrates the wood vessels’ peripheral region extraction comparisons in a ROI with r=30 and different k for wood species Intsia bijuga. The value interval of k is 2,47 with a step length of 5. We can see that a small k will make the wood rays or axial parenchyma remain in vessels’ peripheral region, while a large k will make this region contract.Figure 13
Wood vessels’ peripheral region extraction with differentk values. (a)k=2, (b) k=7, (c) k=12, (d) k=17, (e) k=22, (f) k=27, (g) k=32, (h) k=37, (i) k=42, and (j) k=47.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)Figure14 illustrates the extraction comparisons with k=43 and different r values. The value interval of r is 1,46 with a step length of 5. We can see that a small r will make this region disappear and make wood rays appear, while large r will make this region contract again.Figure 14
Wood vessels’ peripheral region extraction with differentr values. (a) r=1, (b) r=6, (c) r=11, (d) r=16, (e) r=21, (f) r=26, (g) r=31, (h) r=36, (i) r=41, and (j) r=46.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)In summary, we conclude that smallr and k produce more wood rays or axial parenchyma, while large r and k make these tissues disappear while causing vessel peripheral regions to contract. Table 5 gives the best values of r and k for our 6 experimental species, and the corresponding extraction results are illustrated in Figure 15.Table 5
Best values for parametersr and k.
Tree
K
r
Red oak
45
26
Ash wood
30
35
Merbau
43
30
Golden silk teak
45
30
Rose wood
50
45
Limba
20
30Figure 15
Vessel peripheral region extraction with bestr and k for 6 species. (a) Red oak. (b) Ash wood. (c) Merbau. (d) Golden silk teak. (e) Rose wood. (f) Limba.
(a)
(b)
(c)
(d)
(e)
(f)


### 3.5. Incremental Wood Species Identification

To test the incremental identification performance of our scheme, two new wood species were added (i.e., Manchurian Ash (Fraxinus mandshurica Rupr.) and Chinese Oak (Quercus mongolica Fisch.)). It was usually difficult to discriminate these two wood species since they have similar external appearance. The identification results with PCA and three classifiers, BP neural network, SVM, and MD, are illustrated in Figure 16. Again, we can see that the CCR of vessels’ central regions was the lowest in all cases, while the CCR of vessels’ peripheral regions was higher than that of random ROI. In fact, the CCR of vessels’ peripheral regions approaches 90% approximately. Therefore, we conclude that the CCR of vessels’ peripheral regions was stable and best when new wood species samples were added.Figure 16
Spectral classification comparisons for 8 wood species.


### 3.6. Comparisons with Other Works

Finally, the proposed scheme is compared with the two algorithms described in references [15, 16]. Hyperspectral images must be converted to grayscale images before comparison. To this end, two methods can be applied. The first method consists of selecting the waveband of the hyperspectral images using the procedure described in Section 2.5, then performing wavelet-based image fusion on selected images, and finally saving grayscale images after image fusion. This method is detailed below.The retained images after waveband selection were defined asX1,X2,…,Xj,…,Xm,m∈1,2,3,…,10, where the image size of Xj was M×N. Given that each image was derived from a hyperspectral image, image registration of the above images was not necessary. Each image was processed by a wavelet transform based on Haar wavelets. The number of layers for decomposition was 4. The low-frequency wavelet coefficients A1,A2,…,Aj,…,Am and high-frequency wavelet coefficients H1k,H2k,…,Hjk,…,Hmk,k∈1,2,3 of each image were then obtained. The high-frequency wavelet coefficients include horizontal, vertical, and diagonal high-frequency coefficients. The new wavelet coefficients A′ and H′ were then calculated by fusion as(16)A=1mA1+A2+⋯Aj+⋯Am,Hk′=Hjk,s.t.Hjk=maxH1k,H2k,…,Hjk,…,Hmk.These new wavelet coefficients were then used for inverse wavelet transform to obtain the fused digital imageX′. The results of image fusion are presented in Figure 17. Note that the fused image includes the texture features of all of the waveband images.Figure 17
Schematic of image fusion with wavelet transform. (a) BandX1, (b) Band X2, (c) Band X3, and (d) fused image, taking a three-waveband image as an example.
(a)
(b)
(c)
(d)The second aforementioned method consists in applying principal component analysis (PCA) to reduce the dimension of hyperspectral images. The first four principal components are presented in Figure18, but we only retained the first principal component with clear texture.Figure 18
Results of dimensionality reduction from PCA ((a) PC1, (b) PC2, (c) PC3, and (d) PC4).
(a)
(b)
(c)
(d)The method described in reference [15] is a traditional approach to image texture classification applied for the identification of wood species based on GLCM (Gray-Level Co-Occurrence Matrix), which will not be described here. The method described in reference [16] also aims to image texture classification for wood species but based now on IBGLAM (Improved Basic Gray Level Aura Matrix). The corresponding specific steps are described below.According to the descriptions in reference [16], an image needs to be first compressed into a 4-bit grayscale image, which is processed with equation (17). F′x,y refers to an 8-bit grayscale image; floor represents downward rounding; Fx,y refers to a 4-bit grayscale image. The result of grayscale image compression is shown in Figure 19. Note that the image textures are still clear despite the compression of the pixel values of the grayscale image.(17)Fx,y=floorF′x,ymaxF′x,y+1∗16.Figure 19
Result of compression of pixel values ((a) original 8-bit image; (b) compressed 4-bit image).
(a)
(b)Then, the feature matrixB was calculated with the expressions in equation (19). The grayscale image was defined as Fi,j,i∈2,3,…,M−1,j∈2,3,…,N−1, where M and N represent the length and width of the grayscale image, respectively. To make the subscripts of B valid, it is necessary to add 1 to each pixel value of Fi,j before calculation. In this way, we avoid the occurrence of 0, with Fi,j∈1,2,…,16. The matrix B calculated with equation (17) is a 16×16 square matrix.(18)BFi,j,Fi−1,j=BFi,j,Fi−1,j+1,BFi,j,Fi−1,j+1=BFi,j,Fi−1,j+1+1,BFi,j,Fi,j+1=BFi,j,Fi,j+1+1,BFi,j,Fi+1,j+1=BFi,j,Fi+1,j+1+1.To reduce the dimension of the feature vector, the square matrixB was calculated with the following expression:(19)B′=i,j=Bi,j+Bj,i,ifi≠j,Bi,j,ifi=j.Finally, the square matrixB′ was converted according to equation (19), where ⊎ represents the elements in series. The converted vector x was the feature vector with a length of 136. Figure 20 shows the schematic diagram of features. Table 6 shows the processing time for a hyperspectral image when applying the above method.(20)x=⊎B′i,j,i∈1,2,…,16,j∈i,i+1,…,16.Figure 20
IBGLAM feature vector.Table 6
Processing time when applying the methods described in references [15, 16].
Spectral dimension reduction
Texture feature computation
K-L divergence
PCA
IBGLAM
GLCM
42.9767 (s)
0.6414 (s)
0.1698 (s)
0.0329 (s)Figure21 shows the CCR of the vessel peripheral spectral extraction method proposed in this paper in comparison with the IBGLAM and GLCM methods based on support vector machines (SVM) for the identification of eight tree species, including the six tree species in Table 1 and the two tree species in Section 3.5.Figure 21
Comparison of proposed method vs. other methods: (a) vessel peripheral spectral extraction method; (b) IBGLAM and GLCM methods.
(a)
(b)It can be seen from Figure21 that the CCR of the two methods described in references [15, 16] were only around 80% for the samples used in this paper, whereas the CCR of our method proposed in this paper was much higher than that of the other two methods. In addition, when PCA was used for generating grayscale images, the accuracy of image classification with GLCM and IBGLAM was only 70.93% and 73.70%, respectively.


## 3.1. Spectral Difference Analysis for Wood Vessels’ Peripheral Region, Central Region, and Randomly Selected ROI

We may use spectral properties of a specific wood region or tissue to guide wood species identification using both interclass and intraclass differences. It is obvious that large interclass differences and small intraclass differences will produce more reliable results. These two differences can be represented by the spectral curve’s interclass and intraclass distances. Therefore, here, we use the divergence matrix to indicate the spectral curve’s interclass and intraclass differences. We define the intraclass divergence matrix for thejth wood species as follows:(10)Sj=1nj∑k=1njXkj−mjXkj−mjT.Herej=1,2,…,cc=6; Xkj represents the kth sample for the jth species, mj represents the mean value of all samples for the jth species, and nj represents the number of all samples for the jth species. The overall intraclass divergence matrix can be defined as follows:(11)Sw=∑j=1cPjSj,where Pj represents the prior probability for the jth wood species, and it can be estimated by Pj=nj/n, where n represents the overall number of all species samples. Similarly, the interclass divergence matrix is defined as follows:(12)Sb=∑j=1cPjmj−mmj−mT,where m represents the mean value of all species samples. The overall divergence matrix is then defined as follows:(13)St=Sb+Sw.Finally, the pattern divisibility rule is defined as follows:(14)J=trSbtrSw,where trX represents the trace of matrix X. It is obvious that a higher J value will produce a more reliable identification.In terms of a specific wood region or tissue, we considered the wood vessels’ peripheral region, central region, and randomly selected ROI. Spectral extractions in wood vessels’ peripheral and central regions have been discussed in Sections2.4 and 2.5. The size of the randomly selected ROI was 300×300; the mean vector of all spectral vectors in this ROI is the final spectral extraction for a hyperspectral cube sample.The smoothing procedure for the spectral preprocessing and the SNV calibration procedure are usually required before spectral extraction and pattern recognition as they cancel the noise disturbance and waveform offset. This preprocessing process is illustrated in Figure7.Figure 7
Spectral preprocessing procedure: (a) original spectral curve; (b) smoothed spectral curve; (c) spectral curve with SNV processing.
(a)
(b)
(c)Table2 gives the J values for the wood vessels’ peripheral region, central region, and randomly selected ROI in the training samples. It indicates that the spectral divisibility for the wood vessels’ peripheral region was most useful, while that for the wood vessels’ central region was least useful. Since the wood vessels’ central region is hollow, with some resin and probably some filling material, its spectral curve is variable and thus has a large intraclass difference with a small J value.Table 2
J values for different regions of the vessels.
Item
Vessels’ central region
Vessels’ peripheral region
Random ROI
J
1.6506
9.0970
6.6387Finally, the time required was compared for spectral extractions in the wood vessels’ peripheral region, central region, and random ROI. The computer configuration used was as follows: CPU Intel Core I7-6700, internal memory 8G, hard disk 1T, display card AMD Radeon R7200 Series. In these three spectral extractions, one hyperspectral image was cut into 30 ROI of100×100. Time required for processing one ROI is given in Table 3. We can see that the processing speed of random ROI was the fastest while that for the vessels’ central region was the slowest.Table 3
Time requirement comparisons for 3 spectral extractions.
Method
Process
Computing time (s)
Vessels’ peripheral region
G
x
,
y
0.6095
D
x
,
y
0.0056
E
x
,
y
0.0008
Total
0.6199
Vessels’ central region
D
0.7714
K-means cluster
0.2252
Central region selection
0.0087
Total
1.0181
Random ROI
Total
0.0874


## 3.2. Comparisons of Identification Accuracy

Before sending the extracted spectral vectors into classifiers, spectral dimensionality reduction was usually required. In this work, three algorithms, PCA, KPCA, and MDS, were used to reduce dimensionality. In spectral classification, three classifiers, BP neural networks, SVM, and MD, were used for identification comparisons.Firstly, we used the MD classifier to calculate the MD between a test sample and the training sample sets of the 6 hardwood species. The usual nearest neighbor classification was used to judge the test sample’s wood species. The classification results are illustrated in Figure8. We can see that the classification correct ratio (i.e., CCR) was influenced by spectral dimensionality reduction. For instance, if MDS was used, the CCR stabilized only when the spectral dimensionality was larger than 7. On the other hand, if PCA or KPCA was used, the CCR stabilized only when the spectral dimensionality was larger than 3. Moreover, we can see that the CCR was also influenced by the wood region of the spectral collection. In fact, the CCR of the wood vessels’ peripheral region was the highest, approaching 100%, while the CCR of the wood vessels’ central region was the lowest, and its mean value was approximately 70%–80%. The CCR of random ROI fell between the peripheral and central vessel regions and was approximately 90%.Figure 8
Spectral classification comparisons with MD classifier.The classification accuracy was the calculated CCR, which was the overall classification ratio (%) for the 540 test set samples. This classification accuracy is illustrated in Figures8–11 and Table 4.Figure 9
Identification accuracy of BP network for 25 training times (i.e., 25 networks).Figure 10
Spectral classification comparisons with BP neural network classifier.Figure 11
Spectral classification comparisons with SVM classifier.Table 4
Best classification results with different classifiers and spectral dimensionality reductions.
Dimensionality reductions
PCA
MDS
KPCA
Classifiers
Spectral region
CCR (%)
Spectral dimension
CCR (%)
Spectral dimension
CCR (%)
Spectral dimension
BP
Central vessel region
81.78
13
80.42
13
79.4
13
Random ROI
97.80
9
97.57
9
94.25
9
Peripheral vessel region
98.67
7
98.14
8
97.65
9
MD
Central vessel region
81.29
14
78.14
11
66.48
12
Random ROI
88.88
13
91.11
11
94.44
15
Peripheral vessel region
100
8
100
9
98.70
7
SVM
Central vessel region
83.70
15
83.33
14
83.33
13
Random ROI
97.78
7
97.78
8
94.62
12
Peripheral vessel region
99.07
8
99.07
7
98.88
7Secondly, we used the BP neural network classifier. Since the threshold and weight for the network initialization were selected randomly, the classification accuracy of every trained network may fluctuate slightly. Figure9 illustrates the test CCR for the 25 trained networks, while the spectral dimensionality was reduced to 5D by PCA. We observed that the overall CCR was stable, though there were slight differences between each CCR across the 25 networks. Therefore, the mean CCR for 25 networks was calculated and used as the final CCR. Figure 10 illustrates the mean CCR for different spectral dimensionality reductions and different wood regions of spectral collection. We can see that CCRs from vessel central regions were lowest, with mean values of approximately 70%–80%, while those from peripheral regions and random ROI were almost identical, approaching 100%.Thirdly, the SVM classifier was used. The predictive power of SVM depends on the kernel function used. The radial basis kernel function was usually used. When determining the penalty factorC and γ value, the grid search method can be used to find the best parameters. Figure 11 illustrates the identification comparisons with the best parameters C and γ. We can see again that the CCR of central regions was lowest, with a mean of approximately 80%, while those of peripheral regions and ROI were almost identical, approaching 100%. The CCR of random ROI decreased with the application of KPCA.The best identifications with different classifiers and spectral dimensionality reductions are given in Table4.Figures8, 10, and 11 and Table 4 indicated that the CCR of vessels’ central regions was the lowest in all cases and that the CCR of vessels’ peripheral regions was slightly higher than that of random ROI in the case of BP network and SVM. In the case of MD classifier, the CCR of vessels’ peripheral regions was much higher than that of random ROI. Moreover, we obtained good CCR of peripheral regions with a low spectral dimension. For example, the best CCR (100% of vessels’ peripheral region) was achieved using PCA and MD classifier with a spectral dimension of 8 or by using MDS and MD classifier with a spectral dimension of 9.Concerning the spectral dimensionality reductions, good results were obtained using PCA and KPCA with a low spectral dimension. However, if MDS was used, the CCR stabilized only when the spectral dimensionality was larger than 7. In general, PCA was the best choice in terms of identification accuracy and spectral dimension.


## 3.3. Comparisons of Identification Robustness

Wood hyperspectral image acquisition may be influenced by external environmental factors such as temperature, humidity, and illumination. In this section, a simulation experiment on noisy wood hyperspectral image classification was performed. Since the training set of wood images can be acquired in advance in a controlled environment, we only added white noises to the test set. The signal-to-noise ratio (SNR) was used to measure the test set’s image quality.The robustness of our three spectral extraction schemes (i.e. vessels’ central, vessels’ peripheral, random ROI) was compared using an error rate defined in terms of SNR.(15)Error rate=CCRideal−CCRnoisy,where CCRideal and CCRnoisy refer to the classification accuracy for ideal test images and noisy test images, respectively. The falling trend of “error rate” (Figure 12) indicates that the faster the curve falls, the more robust the corresponding spectral extraction scheme. Therefore, we conclude that in the case of BP and MD classifiers, the vessels’ peripheral extraction scheme was the most robust, while for SVM, the vessels’ peripheral scheme was slightly more robust than random ROI and much more robust than the vessels’ central scheme.Figure 12
Different SNR with error rates: (a) SVM; (b) BP; (c) MD.
(a)
(b)
(c)


## 3.4. Discussions of Parameters

In the wood vessels’ peripheral region extraction, the two key parameters are the sizer of structure element A and the threshold k in equation (2). Here we explain the effect of these two parameters and give the best values of r and k for each species.Figure13 illustrates the wood vessels’ peripheral region extraction comparisons in a ROI with r=30 and different k for wood species Intsia bijuga. The value interval of k is 2,47 with a step length of 5. We can see that a small k will make the wood rays or axial parenchyma remain in vessels’ peripheral region, while a large k will make this region contract.Figure 13
Wood vessels’ peripheral region extraction with differentk values. (a)k=2, (b) k=7, (c) k=12, (d) k=17, (e) k=22, (f) k=27, (g) k=32, (h) k=37, (i) k=42, and (j) k=47.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)Figure14 illustrates the extraction comparisons with k=43 and different r values. The value interval of r is 1,46 with a step length of 5. We can see that a small r will make this region disappear and make wood rays appear, while large r will make this region contract again.Figure 14
Wood vessels’ peripheral region extraction with differentr values. (a) r=1, (b) r=6, (c) r=11, (d) r=16, (e) r=21, (f) r=26, (g) r=31, (h) r=36, (i) r=41, and (j) r=46.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)In summary, we conclude that smallr and k produce more wood rays or axial parenchyma, while large r and k make these tissues disappear while causing vessel peripheral regions to contract. Table 5 gives the best values of r and k for our 6 experimental species, and the corresponding extraction results are illustrated in Figure 15.Table 5
Best values for parametersr and k.
Tree
K
r
Red oak
45
26
Ash wood
30
35
Merbau
43
30
Golden silk teak
45
30
Rose wood
50
45
Limba
20
30Figure 15
Vessel peripheral region extraction with bestr and k for 6 species. (a) Red oak. (b) Ash wood. (c) Merbau. (d) Golden silk teak. (e) Rose wood. (f) Limba.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.5. Incremental Wood Species Identification

To test the incremental identification performance of our scheme, two new wood species were added (i.e., Manchurian Ash (Fraxinus mandshurica Rupr.) and Chinese Oak (Quercus mongolica Fisch.)). It was usually difficult to discriminate these two wood species since they have similar external appearance. The identification results with PCA and three classifiers, BP neural network, SVM, and MD, are illustrated in Figure 16. Again, we can see that the CCR of vessels’ central regions was the lowest in all cases, while the CCR of vessels’ peripheral regions was higher than that of random ROI. In fact, the CCR of vessels’ peripheral regions approaches 90% approximately. Therefore, we conclude that the CCR of vessels’ peripheral regions was stable and best when new wood species samples were added.Figure 16
Spectral classification comparisons for 8 wood species.


## 3.6. Comparisons with Other Works

Finally, the proposed scheme is compared with the two algorithms described in references [15, 16]. Hyperspectral images must be converted to grayscale images before comparison. To this end, two methods can be applied. The first method consists of selecting the waveband of the hyperspectral images using the procedure described in Section 2.5, then performing wavelet-based image fusion on selected images, and finally saving grayscale images after image fusion. This method is detailed below.The retained images after waveband selection were defined asX1,X2,…,Xj,…,Xm,m∈1,2,3,…,10, where the image size of Xj was M×N. Given that each image was derived from a hyperspectral image, image registration of the above images was not necessary. Each image was processed by a wavelet transform based on Haar wavelets. The number of layers for decomposition was 4. The low-frequency wavelet coefficients A1,A2,…,Aj,…,Am and high-frequency wavelet coefficients H1k,H2k,…,Hjk,…,Hmk,k∈1,2,3 of each image were then obtained. The high-frequency wavelet coefficients include horizontal, vertical, and diagonal high-frequency coefficients. The new wavelet coefficients A′ and H′ were then calculated by fusion as(16)A=1mA1+A2+⋯Aj+⋯Am,Hk′=Hjk,s.t.Hjk=maxH1k,H2k,…,Hjk,…,Hmk.These new wavelet coefficients were then used for inverse wavelet transform to obtain the fused digital imageX′. The results of image fusion are presented in Figure 17. Note that the fused image includes the texture features of all of the waveband images.Figure 17
Schematic of image fusion with wavelet transform. (a) BandX1, (b) Band X2, (c) Band X3, and (d) fused image, taking a three-waveband image as an example.
(a)
(b)
(c)
(d)The second aforementioned method consists in applying principal component analysis (PCA) to reduce the dimension of hyperspectral images. The first four principal components are presented in Figure18, but we only retained the first principal component with clear texture.Figure 18
Results of dimensionality reduction from PCA ((a) PC1, (b) PC2, (c) PC3, and (d) PC4).
(a)
(b)
(c)
(d)The method described in reference [15] is a traditional approach to image texture classification applied for the identification of wood species based on GLCM (Gray-Level Co-Occurrence Matrix), which will not be described here. The method described in reference [16] also aims to image texture classification for wood species but based now on IBGLAM (Improved Basic Gray Level Aura Matrix). The corresponding specific steps are described below.According to the descriptions in reference [16], an image needs to be first compressed into a 4-bit grayscale image, which is processed with equation (17). F′x,y refers to an 8-bit grayscale image; floor represents downward rounding; Fx,y refers to a 4-bit grayscale image. The result of grayscale image compression is shown in Figure 19. Note that the image textures are still clear despite the compression of the pixel values of the grayscale image.(17)Fx,y=floorF′x,ymaxF′x,y+1∗16.Figure 19
Result of compression of pixel values ((a) original 8-bit image; (b) compressed 4-bit image).
(a)
(b)Then, the feature matrixB was calculated with the expressions in equation (19). The grayscale image was defined as Fi,j,i∈2,3,…,M−1,j∈2,3,…,N−1, where M and N represent the length and width of the grayscale image, respectively. To make the subscripts of B valid, it is necessary to add 1 to each pixel value of Fi,j before calculation. In this way, we avoid the occurrence of 0, with Fi,j∈1,2,…,16. The matrix B calculated with equation (17) is a 16×16 square matrix.(18)BFi,j,Fi−1,j=BFi,j,Fi−1,j+1,BFi,j,Fi−1,j+1=BFi,j,Fi−1,j+1+1,BFi,j,Fi,j+1=BFi,j,Fi,j+1+1,BFi,j,Fi+1,j+1=BFi,j,Fi+1,j+1+1.To reduce the dimension of the feature vector, the square matrixB was calculated with the following expression:(19)B′=i,j=Bi,j+Bj,i,ifi≠j,Bi,j,ifi=j.Finally, the square matrixB′ was converted according to equation (19), where ⊎ represents the elements in series. The converted vector x was the feature vector with a length of 136. Figure 20 shows the schematic diagram of features. Table 6 shows the processing time for a hyperspectral image when applying the above method.(20)x=⊎B′i,j,i∈1,2,…,16,j∈i,i+1,…,16.Figure 20
IBGLAM feature vector.Table 6
Processing time when applying the methods described in references [15, 16].
Spectral dimension reduction
Texture feature computation
K-L divergence
PCA
IBGLAM
GLCM
42.9767 (s)
0.6414 (s)
0.1698 (s)
0.0329 (s)Figure21 shows the CCR of the vessel peripheral spectral extraction method proposed in this paper in comparison with the IBGLAM and GLCM methods based on support vector machines (SVM) for the identification of eight tree species, including the six tree species in Table 1 and the two tree species in Section 3.5.Figure 21
Comparison of proposed method vs. other methods: (a) vessel peripheral spectral extraction method; (b) IBGLAM and GLCM methods.
(a)
(b)It can be seen from Figure21 that the CCR of the two methods described in references [15, 16] were only around 80% for the samples used in this paper, whereas the CCR of our method proposed in this paper was much higher than that of the other two methods. In addition, when PCA was used for generating grayscale images, the accuracy of image classification with GLCM and IBGLAM was only 70.93% and 73.70%, respectively.


## 4. Conclusions

We have proposed a novel hardwood species identification scheme using hardwood microscopic hyperspectral images. The spectral feature of wood vessels is especially considered to investigate the possibility of tree identification to the species level. Specifically, the spectral features of vessel peripheral regions are acquired using morphological operations and those of central regions by use of K-L divergence and an automatic identification scheme. The spectral divisibility rules in Section3.1 and experimental comparisons in Section 3.2 indicate uniformly that the CCR of peripheral regions is best in all cases of spectral dimensionality reduction. Moreover, we can also conclude that wood vessels’ peripheral extraction scheme is robust when wood images are disturbed by white noise.Our work attempted to identify hardwood trees to species level using spectral analysis. Considering the instrument’s cost, we chose the visible/near infrared SOC710VP hyperspectral imager with 128 spectral bands and a spectral range of 372–1038 nm. Spectral data in the visible region may vary due to changes in wood color produced by environmental variation. Therefore, we kept the environment constant during spectral work and wood sample preservation. In addition, 780–900 nm is the near infrared band, which best characterizes different wood properties.The method described in reference [15] is suitable for wood species with large different textures, but the highest CCR was only 75.92% when using the samples provided in this paper. Zamri et al. [16] classified 52 tropical tree species and obtained the highest CCR of 97.01%. However, the main research objects of reference [16] were tropical tree species. These species tend to have large evenly distributed vessels, and most of them are diffuse-porous species that often have clear texture features. The highest CCR was only 82.03% with the method of reference [16] using wood samples provided in this paper. The research objects in this paper include not only tropical tree species, but also subtropical tree species. Moreover, the vessels of these subtropical tree species present different sizes and they are not evenly distributed. Furthermore, the method proposed in this paper does not require rich textures in the cross sections of wood species. As long as the vessels of tree species are clear, the tree species can be classified according to the spectral information of structures around the vessels.



---

*Source: 2039453-2019-06-27.xml*