Abstract
3D shape retrieval is a problem of current interest in different fields, especially in the mechanical engineering domain. According to our knowledge, multifeature based techniques achieve the best performance at present. However, the practicability of those methods is badly limited due to the high computational cost. To improve the retrieval efficiency of 3D CAD model, we propose a novel 3D CAD model retrieval algorithm called VSC_WCO which consists of a new 3D shape descriptor named VSC and Weights Combination Optimization scheme WCO. VSC represents a 3D model with three distance distribution histograms based on vertices classification. The weighted sum of norm distances between corresponding distance histograms of two VSC descriptors is regarded as dissimilarity of two models. For higher retrieval accuracy on a classified 3D model database, WCO is proposed based on Particle Swarm Optimization and existing class information. Experimental results on ESB, PSB, and NTU databases show that the discriminative power of VSC is already comparable to or better than several typical shape descriptors. After WCO is employed, the performance of VSC_WCO is similar to the leading methods by all performance metrics and is much better by computational efficiency.
1. Introduction
The development of computer graphics technology gives a birth to the explosion in the number of 3D models which are now widely used in many diverse applications. Meanwhile, more and more 3D model databases are available on the web, such as Princeton Shape Benchmark (PSB) [1], National Taiwan University Shape Benchmark (NTU) [2], and Engineering Shape Benchmark (ESB) [3]. As Funkhouser et al. [4] predicted, the key question about 3D model has shifted from “How do we construct them?” to “How do we find them?” There is no doubt that searching suitable 3D engineering models from existing resource is extremely beneficial for mechanical design, because () more than 75% of design activity comprises reuse of previous design knowledge to address a new design problem [5] and () modeling a complex 3D model, such as mechanical parts, is still error-prone and time-consuming. Therefore, it is very necessary to develop efficient tools for 3D CAD model retrieval. As a promising technique, the so-called content-based 3D model retrieval, using the characteristics of 3D model itself, has attracted a lot of research interests [6–13]. The content-based approaches are broadly classified into three categories: shape-based, view-based, and multifeature based approaches.
Shape-Based Approaches. This kind of approaches uses the distribution of 3D features to characterize the geometric properties of a 3D model. D2 [14] is a well known shape-based 3D model retrieval method. D2 took samples on the distances between two points on the model surface and then accumulated them into bins of different intervals to form a distance distribution histogram as the model shape descriptor. D2 has many advantages: robustness, rotation invariance, and computational efficiency. However, the discriminative capability of D2 is limited for complex 3D model such as 3D engineering models. Inspired by D2, several new shape-based methods were proposed. Shih et al. [15] investigated a descriptor named Grid D2 (GD2). In GD2, the 3D model is decomposed by a voxel grid. A voxel is regarded as valid if there is a polygonal surface located within this voxel and invalid otherwise. The distribution of distance between two valid voxels instead of two points on the surface is calculated. Different from GD2, the shape distributions in Volume D2 (VD2) [16] are generated from Euclidean distance between randomly selected voxels. Not only surface voxels but also internal ones are involved. In order to reduce the computational complexity of D2, Wu et al. [17] proposed an algorithm called Quick D2 which modifies the sampling phase of D2. Pan et al. [18] proposed a new shape descriptor called Poisson Histogram which can be defined by the following two steps. Firstly, a 3D shape signature was defined by Poisson equation. Secondly, a histogram was constructed by accumulating the values of the defined signature in bins.
View-Based Approaches. The main concept of visual representation in 3D model retrieval is to first convert the 3D model into 2D projection image and then utilize mature image processing techniques to extract various features [19]. The original work was accomplished by Chen et al. [2], who proposed a method named Light Field (LF) descriptor which is a typical view-based approach. The main idea of LF is that if two 3D models are similar, they also look similar from all viewing angles. LF represents a 3D model as a collection of 2D images rendered from uniformly sampled positions on a viewing sphere located around the model. The distance between two descriptors is defined as the minimum L1 distance, taken over all rotations and all pairing of vertices on two dodecahedra. LF had been shown to perform well on the ESB [3] database. In [20], adaptive views clustering (AVC) selected the best characteristic views from more than 320 projected views. The characteristic view selection algorithm is based on an adaptive clustering algorithm and uses statistical model distribution scores to select the optimal number of views. The view-based methods achieved better retrieval performance compared to shape-based ones because they capture more of shape content. However, as pointed out in [2], view-based methods involve more computation compared to the shape-based ones.
Multifeature Based Approaches. According to [3], the shape representations that hold more shape content are better at retrieving more relevant models. Therefore, combining different 3D shape features is an efficient way to enhance the retrieval accuracy [21]. Papadakis et al. [22] presented a novel hybrid 3D model shape descriptor named PANORAMA using a set of panoramic views of a 3D model. A panoramic view is generated by projecting the model to the lateral surface of a cylinder parallel to one of its three principal axes and centered at the centroid of the model. The 3D model is projected to three perpendicular cylinders, each one aligned with one of its principal axes in order to capture the global shape of the model. They used Fourier and wavelet transforms to extract the features for each panoramic view. Leng and Xiong [23] proposed a composite shape descriptor called TUGE which is obtained from a combination of two-view version of depth buffer-based shape descriptor in [24] and GEDT shape descriptor in [4]. Li and Johan [25] proposed a hybrid 3D shape descriptor ZFDR which comprises four features: Zernike moments, Fourier descriptor, depth information and ray-based features.
Multifeatures based methods have more powerful retrieval capability than shape-based and view-based techniques; however, the time consumption is much higher than the latter because of involving more computations. In ZFDR example [25], the average retrieval time on ESB database (contains 866 3D CAD models) is 1.37 sec. If the target database is 100 times bigger than ESB, the average retrieval time will be more than 10 sec. High time cost limits the practicability of multifeature based techniques seriously.
In addition to those methods mentioned above, machine learning technologies are also used to address 3D models retrieval. In [26], 3D models are represented as probability distributions of binary variables on a 3D voxel grid and the proposed method uses a Convolutional Deep Belief Network to learn the distribution of complex 3D shapes across different categories and arbitrary poses from raw CAD data. In [27], the distance between 3D models is computed based on distance histogram features and 3D moment features. Using this distance measure, the relationships between all 3D models in dataset are formulated as a graph structure. A semisupervised learning process is then conducted to estimate the relevance among the 3D models. Leng et al. [28] proposed a 3D model recognition mechanism based on Deep Boltzmann Machines (DBM). The high-level abstraction representation can be obtained from a DBM which is trained by depth images of 3D model and the feature is used in semisupervised classification method. Qin et al. [29] proposed a deep learning approach to automatically classify 3D CAD models according to the mechanical part catalogue. The designed deep neural network classifier is based on the latest machine learning technique, deep learning.
In this paper, to improve the retrieval efficiency of 3D CAD model, we propose a novel retrieval algorithm named VSC_WCO which consists of a new shape-based 3D model descriptor called VSC and Weights Combination Optimization scheme WCO. The VSC descriptor represents a 3D model with three Euclidean distance distribution histograms which are calculated based on vertices classification. The weighted sum of L1 norm distances of corresponding histograms in two VSC descriptors is regarded as the dissimilarity of two models. To avoid risks caused by using fixed weights in VSC distance computation and further improve the retrieval performance on classified 3D model database, WCO is proposed based on Particle Swarm Optimization (PSO) and existing class information of 3D model database. The objective of WCO is searching optimal weights combinations for different query models.
The rest of this paper is organized as follows. The proposed 3D shape descriptor VSC is presented in Section 2. In Section 3, we describe the details of WCO and our 3D CAD models retrieval method VSC_WCO. Experimental results are demonstrated in Section 4. Finally, conclusions and future work are presented in Section 5.
2. 3D Shape Descriptor VSC Based on Vertices Classification
In this section, we define a new shape-based 3D model descriptor named VSC based on classification of vertices. A VSC descriptor consists of three distance distribution histograms which are computed by the following steps: (i) according to Tensor Voting Theory, the vertices of a 3D model are divided into three categories: Face, Sharp Edge, and Corner; (ii) the Euclidean distances between pairs of vertices in the same category are calculated; (iii) the distance distribution histograms of distinct categories are formed.
2.1. Vertices Classification Based on Tensor Voting Theory
Tensor voting theory has great advantages in computer vision tasks such as segmentation and object recognition. According to [30], the normal voting tensor is defined aswhere is vertex of 3D mesh model; is the normal voting tensor of ; denotes a collection of 1-ring neighbor triangles of v; is th 1-ring neighbor triangle in ; represents the unit normal vector of ; is calculated based on [30]:where area () is area of triangle ; is maximum area among ; denotes barycenter of triangle ; represents the position of .
Because is a rank-3 positive semidefinite matrix, it can be diagonalized as follows: where , , and are the corresponding unit eigenvectors of , , and (), respectively. According to the eigenvalues, vertices can be divided into Face, Sharp Edge, and Corner as follows [31]:(i)Face: is dominant; and are close to 0.(ii)Sharp Edge: and are dominant; is close to 0.(iii)Corner: , , and are approximately equal.
Figure 1 shows an example of vertices classification; (a) is a 3D CAD model called advgr01 in ESB. Marked points in (b), (c), and (d) denote the vertices of Face, Sharp Edge, and Corner, respectively. The numbers of vertices in the distinct categories are 670, 127, and 16, respectively.
(a)
(b)
(c)
(d)
2.2. Distance Distribution Histogram
After classification of vertices, the Euclidean distances between any two vertices in the distinct categories are measured. To eliminate the influence of quantity of points, the distances histogram containing n bins is defined as where is the number of distances within the range of the th bin; is the quantity of each type of vertices; denotes a histogram constructed by counting how many distances fall into each bin. The width of bin is determined bywhere and are maximum and minimum distance between pairs of vertices in the same class; n is the number of bins.
Three distance distribution histograms compose the 3D shape descriptor VSC which is defined as follows: where , , and denote distance distribution histograms of vertices in Face, Sharp Edge, and Corner, respectively. Table 1 is a comparison between feature extraction results of D2 and VSC for six models of ESB database, the first two models are taken from different classes of ESB; the third and fourth models are taken from ESB∖Flat-Thin Wall Components∖Back Doors; and the fifth and sixth are taken from ESB∖Solid Of Revolution∖Gear-like Parts. As shown in Table 1, the first and second models are completely different in shape, but their D2 descriptors are very similar. Different to D2, those two models are quite different in the view of VSC. The third and fourth models are taken from the same class and their shapes are very similar. There are almost no differences between the VSC descriptors of those two models, so do the D2 descriptors. The fifth and sixth models are taken from the same class also, but there is a little difference between their shapes. Through comparing two VSC descriptors of those two models, it is clear that the VSC descriptor reflects differences between two models successfully. Generally speaking, due to the subdivision of distances between sample points, the shape discrimination capacity of VSC is better than D2.
2.3. Distance between Two VSC Descriptors
After VSC descriptors are constructed, the similarity comparison between two 3D models is mapped into the comparison of VSC descriptors. The distance between two VSC descriptors is defined as the weighted sum of L1 norm distances of corresponding histograms:where and denote two 3D models; and are VSC descriptors of and , respectively; L1 is the L1 norm distance measure; , , and are weights. The three distance histograms reflect the different shape characteristics of a model and they have the same contribution for VSC distance computation. Therefore, we linearly combine them. In addition, the distances between the corresponding histograms fall in the same range; as such, we assign the same weight for each histogram: . For two 3D models, the distance between corresponding VSC descriptor is regarded as the dissimilarity, and the smaller the VSC distance is, the more similar they are. According to (7), the VSC distance between first two models shown in Table 1 is 1.71.
For a query model , different weights combinations used in (7) will produce different dissimilarities between and models in database and lead to different retrieval results eventually. The influence of weights combination against retrieval result provides an opportunity to achieve more desirable retrieval result on a classified 3D model database (details are discussed in the next section).
3. 3D CAD Model Retrieval Algorithm VSC_WCO Based on VSC and Weights Combination Optimization Scheme WCO
In this section, to further improve retrieval performance in a classified 3D model database, we propose a 3D CAD model retrieval algorithm called VSC_WCO based on 3D shape descriptor VSC presented in Section 2 and Weights Combination Optimization scheme WCO. In order to avoid the risks caused by using fixed weights combination in VSC distance computation, WCO takes into account the class information of 3D model database and utilizes Particle Swarm Optimization (PSO) to search optimal weights combination for each class.
3.1. Weights Combination Optimization
In order to avoid the risks caused by fixed weights and improve retrieval performance, we define Weights Combination Optimization method named WCO which takes into account the class information of 3D model database and utilizes Particle Swarm Optimization (PSO) [32] to search the optimal weights combination for each class. The motivations of WCO are described as follows:(1)For a query model q, different weights combinations used in (7) will produce different dissimilarities between q and models in database and lead to different retrieval results eventually. The influence of weights combination against retrieval result provides an opportunity to achieve more desirable retrieval result by searching more suitable weights combination for q.(2)In fact, those three distance histograms , , and in the proposed shape descriptor VSC depict a 3D model from different aspects, respectively, and they can be regarded as three single 3D shape descriptors. It is demonstrated in [3] that no single 3D model descriptor performs well on all models and different descriptors have different strengths and weaknesses. According to this conclusion, fixed weights combination used in (7) for VSC distance computation is not suitable for all models.
In this subsection, we introduce the PSO first and then describe the details of WCO.
3.1.1. Particle Swarm Optimization (PSO)
PSO has several advantages over other artificial intelligence algorithms. For example, it is better at global optimization and is easier to apply to multiple-objective problems [28]. PSO is initialized with a population of random potential solutions and the algorithm searches the optimal solution according to its performance. The goodness of a particle is determined by a function called fitness function in PSO. The fitness function takes position of a particle as input and returns a single number which denotes the goodness of the particle for the optimization problem.
Consider a group of N particles that are searching a global optimal solution in a D dimensions space. The position and velocity of (th particle) can be expressed aswhere and denote the velocity of in iterations and ; and represent the positions of in iterations and ; is the best position of until iteration ; is the best position among all the particles until iteration ; and are the personal learning coefficient and the social learning factor, respectively; and are random numbers in the range of (0, 1); is an inertia factor from 0.8 to 1.2.
In WCO, PSO is used for searching the optimal weights combination for each class in a classified 3D model database; thus the position and velocity of each particle consist of three numbers, respectively. In order to avoid negative numbers during the optimization process, the relationship between positions and weights is expressed as where is the th weight and is the th component of particle’s position.
It is obvious that if the weights combination used in (7) is more reasonable to query model , the VSC distance between q and all models in database will be more accurate and the retrieval performance will be better. In other words, the more reasonable weights combination, the better retrieval performance. Thus, the retrieval performance achieved by a particle can be used as its fitness. We employ the performance metrics First Tier (FT) [1] to evaluate the retrieval performance of each particle. FT is expressed aswhere R represents a retrieval result which is formed by sorting all models in ascending order based on their VSC distances from query model ; denotes the class of and is the cardinality of ; is the number of models of in top list of .
3.1.2. Implementation Process of WCO
By controlling weight alterations, one ensures the evaluation of retrieval result better, which is the main idea of WCO. The optimization process of WCO is summarized in Algorithm 1 and notations used in Algorithm 1 are explained in the Notations. According to Algorithm 1, the time and space complexity of WCO are and .
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
3.2. Implementation of VSC_WCO
To improve the retrieval performance on a classified 3D CAD model database, we propose a new 3D CAD model retrieval algorithm named VSC_WCO based on VSC descriptor described in Section 2 and Weights Combination Optimization scheme WCO presented in Section 3.1. Our VSC_WCO is composed of two parts: online and offline parts, which are described as follows and the whole process of VSC_WCO is illustrated in Figure 2.
Offline(1)VSC Feature Extraction. We extract the VSC shape descriptors of all the models in Train and Test sets, based on the method in Section 2.(2)Optimal Weights Calculation in Train Set. We calculate the optimal weights combinations (OWC) for all classes of Train set, based on WCO described in Section 3.1. All optimal weights combinations are stored in a database with the corresponding class information.
Online(1)VSC Feature Extraction. We extract the VSC shape descriptor of query model q, based on the method in Section 2.(2)Determine the Most Similar Class (MSC) of q in Train Set. First, we compute the VSC distance between q and models in Train set according to (7) using weights combination: , , and , and then we select the nearest model’s class as the most similar class of q, which is denoted as MSCq.(3)Determine the Optimal Weights Combination (OWC) of q according to MSCq. The optimal weights combination corresponding to q’s MSC is selected from database, which is denoted as OWCq.(4)VSC Distances Computation in Test Set. We compute the VSC distances between q and models in Test set according to (7) using weights combination OWCq.(5)Ranking and Output. Sort all the models in Test set in ascending order based on their VSC distances computed in step () and output retrieval lists accordingly.
4. Experiments
In this paper, all retrieval experiments are implemented in MATLAB R2011b and performed in a PC with configuration: CPU: Intel Pentium Dual-Core [email protected] GHz; memory: 2.0 GB; OS: Windows XP. To investigate the performance of VSC_WCO for 3D engineering and generic models, we selected the following three representative standard benchmark databases:(i)Engineer Shape Benchmark (ESB) [3]. It is developed by Purdue University for evaluating the search methods relevant to the mechanical engineering domain. There are 866 3D engineering models in ESB and those models are classified into three superclasses, namely, Flat-Thin Wall Component, Rectangular-Cubic, and Solids of Revolution. Within each superclass, models are further classified into clusters of similar shapes. Figure 3 shows some example views of 3D models in ESB. In order to ensure the usefulness of weights combinations calculated by WCO, we equally divide ESB database into two parts: Train and Test sets, the former is used for WCO and the latter is used as target database for retrieval experiments.(ii)Princeton Shape Benchmark (PSB) [1]. It contains 1814 3D models totally, which are classified into Test and Train parts. In our experiments, Train database is used for WCO and Test database is used for retrieval experiments.(iii)National Taiwan University (NTU) Database [2]. NTU database contains 1,833 3D models and only 549 3D models are grouped into 47 classes and the remaining 1,284 models are assigned as the “miscellaneous.” Thus, we only use the 549 classified ones for our experiments.
4.1. Retrieval Performance Evaluation Metrics
To comprehensively evaluate the 3D model retrieval results, we employ five metrics including precision-recall curve, Nearest Neighbor (NN), First Tier (FT), Second Tier (ST), and Discounted Cumulative Gain (DCG) [1]. Precision indicates how much percentage of the top K models belongs to the same class as the query model while recall means how much percentage of a class has been retrieved among the top K retrieval list. The precision-recall curve comprehensively demonstrates retrieval performance, which is assessed in terms of average recall and average precision. NN is defined as the percentage of the closest matches that are relevant models. FT is the recall of the top list, where is the cardinality of the relevant class of the query model. ST is the recall of the top 2() list. DCG is defined as the summed weighted value related to it which combines precision and recall as well as ranking positions.
4.2. Comparison Algorithms
To compare the performance of our retrieval algorithm VSC_WCO, we consider the following four algorithms:(i)D2: a classic shape-based retrieval method and the most related method of VSC. We implemented D2 based on the original paper [14].(ii)LF: a typical view-based retrieval method, which performed well on ESB [3]. We performed experiments based on their provided execution files [2].(iii)ZFDR: a state-of-the-art hybrid retrieval algorithm, which integrates a 3D model’s both visual and geometric information from different aspects. For the retrieval performance of ZFDR on ESB, PSB, and NTU, we refer to [25].(iv)PANORAMA: a state-of-the-art hybrid retrieval algorithm, which utilizes a set of panoramic views. According to Li and Johan [25], PANORAMA is arguably the best retrieval method reported to date. For the retrieval performance of PANORAMA on ESB, PSB, and NTU, we refer to [22, 25].
4.3. Experiments Results and Analysis
In our experiments, every model in Test set is used as query model and to avoid bias we exclude this model from Test set when computing the VSC distance. Figures 4 and 5 and Table 2 compare the performance of our VSC_WCO algorithm and the above-mentioned four methods. To evaluate the effectiveness of WCO, for comparison, we also list the performances when using only the VSC shape descriptor. As can be seen in Figures 4 and 5 and Table 2, firstly, our shape descriptor VSC is better than D2 and comparable to LF with lower computational complexity. Secondly, after applying the WCO, VSC_WCO achieves significantly better performance than VSC. On ESB, VSC_WCO exceeds VSC in NN, FT, ST, and DCG by 2.4%, 13.4%, 12.2%, and 6.1%. On PSB, those metrics were increased by 2.9%, 16.7%, 14.8%, and 16.1%. The corresponding increments are 3.1%, 16.2%, 17.4%, and 13.1% on NTU. This increase indicates that fixed weights combination for VSC distance computation is not suitable for all query models and the optimal weights combinations produced by WCO are useful to improve the retrieval performance of VSC. Finally, the performance of VSC_WCO is substantially the same as that achieved by the state-of-the-art algorithms, such as PANORAMA and ZFDR. Furthermore, the average retrieval time of VSC_WCO on ESB, PSB, and NTU database are less than the latter. Table 3 lists the timings of VSC_WCO on ESB, PSB, and NTU databases. In Table 3, denotes the time consumption of VSC extraction of query model; denotes the time consumption of finding the most similar class of query model in Train set; is time consumption of searching similar models in Test set. On a common configuration computer (CPU: Intel Pentium Dual-Core [email protected] GHz; memory: 2.0 GB), the average retrieval time of VSC_WCO is less than 0.5 s which enables real-time retrieval from large repositories. This is a significant advantage of the VSC_WCO as it can be rendered more efficient by reducing its storage requirements and time complexity. Altogether, VSC_WCO has considerable retrieval performance and is very efficient in computation and comparison time.
In Table 4, we provide a few examples of queries and the retrieved results to first 8 retrieval models from the ESB database using D2, VSC, and VSC_WCO. As illustrated in Table 4, the performance of VSC_WCO is better than VSC and D2 in terms of both quantity and rankings of relevant models.
4.4. Limitations
Our 3D CAD model retrieval algorithm VSC_WCO has achieved good performance on ESB database. However, there are still some limitations. Firstly, since the class information of 3D model database is a key condition for the proposed Weights Combination Optimization method WCO, VSC_WCO only can be used to the already classified 3D model database. Secondly, the optimal weights combination to a query model q is determined by the most similar class of q; therefore, if the most similar class is not the target class, then our method VSC_WCO will fail to get the right models.
5. Conclusions and Future Work
In this paper, to improve the retrieval efficiency for 3D CAD models in a classified 3D model database, we have proposed a 3D model retrieval method called VSC_WCO which is based on the proposed 3D shape descriptor VSC and Weights Combination Optimization scheme WCO. VSC descriptor represents a 3D model with three distance histograms which is computed by vertices classification based on Tensor Theory. WCO utilizes the existing class information and Particle Swarm Optimization (PSO) to search optimal weights combinations for all classes. The retrieval experiments are conducted on Engineering Shape Benchmark (ESB). We equally divide the ESB database into two parts, one is used as Train set for WCO to compute the optimal weights combinations for every class in ESB, and the other one is used as Test set for retrieval experiment. Several metrics are selected to evaluate the retrieval capacity of our algorithm. Retrieval results on ESB demonstrated that () the retrieval performance of the proposed 3D shape descriptor VSC is much better than shape-based descriptor D2 and comparable to view-based descriptor Light Field with lower computational complexity and, () after WCO is employed, the performance of VSC_WCO is close to the state-of-the-art algorithms, such as PANORAMA and ZFDR. Furthermore, the average retrieval time of VSC_WCO on ESB database is much lower than the latter.
Further research will focus on the following aspects: () in order to improve the accuracy of the most similar class of query model, we need to enhance the discrimination capacity of proposed shape descriptor VSC; () develop a 3D model clustering algorithm for applying our method VSC_WCO on unclassified 3D model databases.
Notations
| : | Classified 3D models database which has several classes |
| : | Number of models in |
| : | Number of classes in |
| : | The th class of |
| : | Number of models in |
| : | The th model of |
| : | Number of particles in PSO |
| : | The th particle of PSO |
| : | Query model |
| : | VSC distance between and |
| Ascending : | Sort all the models in in ascending order based on VSC distances |
| : | Retrieval result when is used as query model |
| : | Maximum number of iterations |
| : | Current number of iterations |
| : | Fitness of in the th iteration |
| : | Retrieval performance of when is used as query model, measured by metrics First Tier which is expressed as (13) |
| : | The best position achieved by so far |
| : | Three components of ’s position |
| : | Fitness of |
| : | The best position achieved by any particle so far |
| : | Three components of best’s position |
| : | Fitness of best. |
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This study was supported by Natural Science Foundation of China (51001121), Fundamental Research Funds for the Center Universities and Program for Innovation Research of Science in Harbin Institute of Technology (PIRS OF HIT Q201503), China Postdoctoral Science Foundation funded project (2015M581450), and Program for Innovation Research of Harbin City (2015RAQXJ090).