The two-dimensional transfer functions (TFs) designed based on intensity-gradient magnitude (IGM) histogram are effective tools for the visualization and exploration of 3D volume data. However, traditional design methods usually depend on multiple times of trial-and-error. We propose a novel method for the automatic generation of transfer functions by performing the affinity propagation (AP) clustering algorithm on the IGM histogram. Compared with previous clustering algorithms that were employed in volume visualization, the AP clustering algorithm has much faster convergence speed and can achieve more accurate clustering results. In order to obtain meaningful clustering results, we introduce two similarity measurements: IGM similarity and spatial similarity. These two similarity measurements can effectively bring the voxels of the same tissue together and differentiate the voxels of different tissues so that the generated TFs can assign different optical properties to different tissues. Before performing the clustering algorithm on the IGM histogram, we propose to remove noisy voxels based on the spatial information of voxels. Our method does not require users to input the number of clusters, and the classification and visualization process is automatic and efficient. Experiments on various datasets demonstrate the effectiveness of the proposed method.

1. Introduction

Direct volume rendering [1] has been widely used in many fields, particularly for the visualization of medical data with a variety of modalities such as CT, MRI, and ultrasound. It projects three-dimensional (3D) volumetric data to a two-dimensional (2D) screen to facilitate observation and exploration. By using appropriate transfer functions (TFs), which map the voxel properties (e.g., gray scale and gradient magnitude) to optical properties (e.g., transparency and color), the major structures of volume data can be revealed. Designing effective TFs is a must for useful visualization of volumetric medical data, especially for clinical diagnosis and treatment. However, it remains a challenging task for radiologists and physicians, as it usually requires them to acquire technical knowledge on rendering techniques. Unfortunately, doctors do not have enough time to acquire the related knowledge and skills. In addition, the complicated interactions in traditional direct volume rendering systems prohibit their application in clinical practice. In this regard, developing automatic methods for TFs generation is important for medical data visualization.

TFs typically operate on the value domain of input volume and its derived feature domains (such as gradient) to separate different structures. To assist users’ exploration, the feature domains are always decomposed into several meaningful clusters and users can explore the data based on these clusters. The voxels’ scale values and their gradients are the most commonly used feature domains. Kindlmann and Durkin [2] projected the scale values and gradients into a 2D space to make up the intensity-gradient magnitude (IGM) histogram and then designed TFs based on it by using various manipulation widgets (such as rectangles, triangles, and ellipses). Later, IGM histogram was often used in direct volume rendering and has been implemented in several popular visualization packages, such as VisIt [3], Voreen [4], and ImageVis3D [5]. However, most implementations involved a lot of manual operations to achieve satisfactory results, which made them not applicable for doctors without special knowledge of IGM histogram. In recent years, Wang et al. [6] proposed a volume exploration with the Gaussian mixture model and generated a set of suggestive elliptical transfer functions semiautomatically. Nonetheless, this method is not effective when the dataset is complicated, as the widgets are too regular to describe multidimensional features with complex shapes. Ip et al. [7] presented a semiautomatic approach to detect embedded features and spatial structures by visually segmenting the IGM histogram of a volumetric dataset. However, this method required too many operations to annul undesired structures and achieve satisfactory results and hence was still not suitable for radiologists and physicians.

To overcome these shortcomings, we propose a novel automatic TFs design method based on IGM histogram. Our method is based on the observation that voxels belonging to the same structure usually have similar intensity and gradient magnitude and are located together in the volume. In this case, we employ affinity propagation clustering algorithm to classify the scattering points in the IGM histogram into clusters by defining two novel similarity measurements: IGM similarity and spatial similarity. These two similarity measurements can effectively bring the voxels of the same tissue together and differentiate the voxels of different tissues so that the generated TFs can assign different optical properties to different tissues, making the visualization results unambiguous and easy to be used in diagnosis and treatment. Note that as we do not need to assign a cluster number to the affinity propagation clustering algorithm in advance, the clustering process is fully automatic. In addition, before performing clustering, we eliminate noisy voxels by leveraging the spatial information of the input volume data. We also provide a simple interaction interface to allow users further polish the clustering results. The TFs are automatically generated based on the clustering results. Experiments on various datasets demonstrate the effectiveness of the proposed method.

The organization of this paper is as follows. We presents related work in Section 2. Section 3 describes our method in detail. Section 4 provides experimental results. Conclusions are drawn in Section 6.

Although direct volume rendering technique has been widely used, transfer function design remains a challenging task. The design of TFs is often classified into two categories: image-centric and data-centric [8], depending on whether they derive their parameters from the resulting images or original data.

Image-centric methods typically adjusted TFs by searching for the optimal rendered image in the rendering space. Thus, the TF’s generation is considered as a process of parameter optimization. In early investigations, genetic algorithms were employed to select the most satisfactory rendering result [9]. However, depth information cannot be sufficiently considered in these algorithms. To enhance the depth information of rendering images, Chan et al. [10] proposed a novel metric based on psychological principles, and Zheng et al. [11] proposed an effective design by optimizing an energy function based on quantitative perception models. Although image-centric methods are intuitive and easy to implement, the final rendering results generated by these methods heavily depend on the initial result. In addition, a lot of parameters need to be tuned in these methods, which is unacceptable for medical data visualization.

In contrast with image-centric methods, data-centric TFs define visual properties based on the information of voxels. In early studies, TFs were designed based only on the histogram of voxels’ scalar values (these TFs are usually called 1D TFs). However, 1D TFs have inherent difficulties in identifying different materials with similar intensities. Then 2D TFs based on intensity values and gradient magnitudes were proposed, which were more effective in detecting multiple materials and their boundaries [2]. Besides using derivative properties of scalar value, other effective metrics were also incorporated into the feature space, such as curvature [12], feature size [13], and ambient occlusion [14]. In order to highlight important structures of volume, some methods based on visibility were introduced [15] and information divergence between visibility distribution and target distribution was optimized to generate automatic TFs [1618]. Recently, Cai et al. have proposed a rule-based transfer function based on the local frequency distribution of a voxels neighborhood [19]. The structures of user’s interest in a volume dataset can be characterized by using the various feature space. However, higher dimensions of the feature space add more complication in transfer function design.

Recently, several semiautomatic transfer function designs have been proposed to reduce the number of degrees of freedom in transfer function design. Tzeng and Ma proposed a method based on material classes extracted from the cluster space using the ISODATA technique to generate TFs [20]. Some other machine learning based methods have also been employed to analyze the feature spaces. For example, MacIejewski et al. proposed a nonparametric kernel density estimation method onto the IGM feature space in order to extract feature patterns [21]. Furthermore, Wang et al. applied the Morse theory to automatically decompose the IGM feature space into a set of valley cells [22]. In order to display the surface information of the volume data better, Šereda et al. proposed LH histogram which uses the lower FL and higher FH intensity values of the two materials that form a boundary as the two axes to design the transfer function [23]. Based on LH histogram, Nguyen et al. proposed to oversegment LH histogram using mean shift clustering and then group similar voxels by using hierarchical clustering [24]. Roettger et al. incorporated spatial information into the IGM histogram and performed clustering based on the barycenter and variance of each bin [25]. Ip et al. introduced a hierarchical volume exploration scheme based on normalized-cut algorithm [7].

The affinity propagation (AP) clustering algorithm is an effective clustering algorithm proposed by Frey and Dueck [26]. It clusters based on similarity of data and has no requirement for the symmetry of the similarity matrix of data, giving it a relatively wide application. AP algorithm regards all data points as potential cluster centers simultaneously and produces high-quality cluster centers by using message transfer between the data points. Compared with -means clustering algorithms [27], AP algorithms do not need to prespecify the number of clusters . Although the hierarchal clustering algorithm [28] also does not require the number of clusters, its classification results [23, 24] are sensitive to the initial clustering result while AP algorithm clustering resolves this issue by improving the reliability of clustering result. As mentioned earlier, the normalized cut algorithm [29] was also employed in volume visualization [7, 17], which overcomes the disadvantages of the -means algorithm for its ability to identify nonconvex distribution clustering, but it requires the solution of matrix eigenvalues and eigenvectors, which is quite computation-intensive and time-consuming. In contrast, convergence speed of AP algorithm is much faster, and its clustering results are more accurate [26].

In this paper, we propose an automatic method to classify the volume data not only considering the intensity value and gradient magnitude but also taking the spatial information of voxels into consideration. Then the TFs are automatically generated based on the clustering results.

3. Method

Figure 1 shows the workflow of the proposed method, which is composed of three steps. First, we calculate the gradient at each voxel. Together with the intensity, the IGM histogram of the volume is constructed. We call a set of voxels with the same intensity and gradient value a bin. We compute the mean and variance of position of the voxels within each bin. The bin with a larger variance of position than a specified threshold is eliminated from the histogram, as it can be considered as noise or nonsignificant tissues located at the boundaries of the volume. Second, we cluster the remainder bins in the histogram using affinity propagation algorithm according to both the intensity and gradient information and the spatial information of voxels. In addition, we provide a simple interaction interface to allow users to further polish the clustering results in order to achieve more desired visualization. Finally, the transfer functions are automatically generated based on the clustering results, and the visualization of the volume can be obtained from the transfer function.

3.1. IGM Histogram

The IGM histogram is a useful tool for exploring the volume dataset [2, 7, 25]. It is generally represented as a 2D scatter plot. To figure out the IGM histogram from a given 3D intensity field, we first compute the gradient magnitude at each voxel. Then, the intensity and gradient magnitude is set as the horizontal and vertical axis of the scatter plot. Each point in the scatter plot represents a bin of voxels with the same intensity and gradient value. The brightness value of each point is associated with the number of the corresponding voxels. To avoid significant changes of the magnitude of the brightness, we set it as a logarithmic function of the number of voxels. If all bins are fed into the following clustering algorithm, the generated transfer function obtained from the clustering results as well as the final visualization may be greatly affected by a great deal of noise. In this regard, we should try our best to maintain the bins of important tissues while eliminating nonsignificant and noisy bins as much as possible to facilitate the histogram clustering. It is observed that the spatial location of important tissues and organs are concentrated. Based on this observation, we employ a simple and efficient scheme to refine the IGM histogram considering the spatial information of voxels. We first calculate the mean position of voxels and the variance of voxel positions bywhere is the number of voxels in the bin, is the voxel set, and is the position of voxel . If is larger than a threshold value, we remove the bin from the histogram. Figure 2 shows the intensity-gradient histogram of abdomen dataset (see Figure 6) with different threshold values, where gray means the removed bins based on the prespecified threshold while the red means the left bins for clustering. More results can be found in Section 4.

3.2. Clustering of IGM Histogram
3.2.1. Basics of Affinity Propagation Clustering

For the completeness of this paper, we introduce the basics of the affinity propagation (AP) clustering algorithm here and readers can refer to [26] for more details. AP clustering algorithm initially regards all the data points as potential cluster centers. Messages between data points then begin to transmit iteratively to maximize a fitness function until an optimal set of exemplars and relative clusters emerge. The similarity between data point and data point , notated as , is generally calculated by negative Euclidean distance between them. The value of similarity increases when the distance between the two points decreases. The elements on the diagonal of the similarity matrix are defined as “preference.” Initially, the value of is set as the median of the similarities between data point and all other data points.

The most important messages in AP algorithms are responsibility and availability. The responsibility is the message sent from data point to candidate cluster center , which reflects whether is suitable to be the cluster center of data point , taking into account other potential exemplars for . The availability is the message sent from candidate cluster center to data point , which reflects whether data point should choose as its cluster center. The alternating renewal process of these two kinds of messages is the core of the AP algorithm.

The initial values of responsibility and availability are set to be . For data point , the responsibility and the availability are updated in an iterative way asThe self-availability is updated in a slightly different way asUpon convergence, the exemplar for each data point will be chosen by maximizing the following criterion:Finally, we obtained a set of exemplars and corresponding clusters. The similarity matrix plays an important role in AP algorithm, as the calculation of both the responsibility and availability depends on the similarity matrix. Traditionally, the similarity matrix in AP algorithm is calculated based on the negative value of the Euclidean distance between two data points. However, it is not sufficient for our application, where the data points are bins on the IGM histogram and their Euclidean distance only contain the intensity and gradient magnitude information. As mentioned, the spatial information of voxels are also quite important to achieve desired visualization results. To this end, we design our own similarity measurement which integrates both intensity-gradient information and spatial information of voxels.

3.2.2. Similarity Measurement

The similarity measurement employed in our application should be designed to bring the voxels of the same tissues together and differentiate the voxels of different tissues so that the generated transfer function can assign different opacity and color to different tissues, making the visualization results unambiguous and easy to be used in diagnosis. As mentioned, statistically, the voxels of the same tissue have similar intensity and gradient magnitude and are located closely within the volume. However, most previous works only use the information of IGM histogram and ignore the location information of voxels, which may influence the visualization result. In order to take full advantage of both IGM and spatial information, we integrate IGM similarity measurement and spatial similarity measurement into the AP clustering algorithm to achieve a more effective clustering results for more appealing visualization of the volumetric data.

IGM Similarity. Usually, close bins in the IGM histogram have similar intensity and gradient magnitude. In this regard, we employ Euclidean distance to measure the IGM similarity of two bins and :where and are the intensity magnitude of the two bins and and are the gradient magnitude of the two bins. We normalize and obtain the IGM similarity:where and are the maximum and minimum value of .

Spatial Similarity. The second similarity measurement is designed to group bins that are neighbored in the volume. We call it spatial similarity. We evaluate the spatial similarity between two bins using the number of direct neighborhood relations between the two bins; we employed the method introduced in [23] to calculate the number. Given two bins and , the number of direct neighborhood relations between and can be calculated bywhere represents voxels in ; represents voxels in ; and is a Boolean function and we set it as if voxel is a neighbor of voxel and otherwise. Note that we define .

The total number of neighborhood relations of bin iswhere represents all other bins in the IGM histogram. Then the spatial similarity between two bins () is computed by taking the maximum of the normalized sum of neighborhood relations:where and means that there no any neighborhood voxels of .

By integrating the IGM similarity and the spatial similarity, the similarity measurement employed in the AP clustering algorithm can be defined aswhere the constant and constant are the weights of the two similarity measurements. We set and in our experiments.

3.2.3. Affinity Propagation Clustering on IGM Histogram

Once the similarity measurement is defined, we can figure out the similarity matrix using (11) and then feed it into the AP clustering algorithm. In order to avoid numerical oscillations, we introduce a damping parameter into (3)–(4), and these equations can be reformulated aswhere the superscripts and represent the number of iterations. The damping parameter is between and , and we set it as a default value in our experiments. The message-passing procedure may be terminated after changes in the messages fall below a threshold.

The whole AP clustering can be summarized as follows:(1)Calculate the gradient magnitude of each voxel.(2)Calculate the variance of each bin.(3)Remove the noisy bins according to the threshold.(4)Calculate the IGM similarity and spatial similarity and obtain the similarity matrix according to (11).(5)The availability and the responsibility matrixes are initialized as .(6)Do the following until changes in the messages fall below a threshold.(i)Calculate the availability and responsibility .(ii)Update and according to (12) and (13).(7)Output the AP clustering result.

3.3. Generation of Transfer Function
3.3.1. Opacity Transfer Function

After performing AP clustering on IGM histogram, the voxels in the volume data are classified into a set of clusters. We first assign an opacity value to each cluster based on the average distance between the voxels of the cluster and the center point of the volume. In order to achieve more appealing visualization, we further refine the opacity value of each voxel in the cluster according to the voxel’s gradient value, which usually indicates whether the voxel is near the boundary of the cluster or not.

A cluster is more likely to occlude other clusters when it is located closed to the boundary of the volume and farther to the center point of it. The average distance between the voxels of the cluster and the center point of the volume can be calculated bywhere is the volume’s center; is the number of voxels in cluster ; and represents voxels of the cluster . We employ the average distance to define the opacity value of the cluster:where the [, ] is predefined opacity range and and are the maximum and minimum average distance between clusters and volume’s center.

For clear-cut visualization, the boundaries of a cluster are much more important than the internal region of the cluster, as boundaries determine the shape of a cluster (i.e., an organ or tissue). In addition, the depiction of boundaries can also enhance the perception of local occlusions, which are quite essential in diagnosis and surgical planning. In this regard, we should make the boundaries of a cluster visually distinct by differentiating the voxels closer to the boundaries and the voxels located in the internal region. In order to achieve this effect, we refine the opacity value at each voxel according to the gradient magnitude of this voxel:where is the opacity value of the voxel; is the opacity value of the cluster ; is the gradient magnitude of the voxel; is the maximum gradient magnitude in the cluster ; and is a parameter that determines how to map the gradient value to the opacity value (we set as or in our experiments).

3.3.2. Color Transfer Function

Color is also an important tool for distinguishing different organs or structures. We assign voxels in a cluster the same color. In our implementation, the HSV model was applied for color mapping as it is a perception-enhanced color model. We divided the space of hue () equally according to the number of the clusters first. Then, the saturation () and value () are set to constants as and . Specifically, for a voxel which is classified into the cluster , the hue can be computed aswhere is the number of clusters.

3.3.3. Interaction Interface

For flexibility, we also provide an interaction interface for our visualization system. In traditional visualization systems that employ IGM histogram to design transfer functions, users are usually requested to select voxels for visualization in the IGM histogram by manipulating various widgets, such as rectangle, triangle, and ellipse. In this case, users are often required to have a good understanding of the IGM histogram so that the widgets can be put in appropriate places for meaningful visualization. It is not an easy task for doctors. In contrast, in our visualization system, as we segment the IGM histogram into several clusters, users can interact with the clusters instead of scattered points in the IGM histogram, which makes the interaction more intuitive and effective. For example, users can merge two or more clusters by setting the clusters with the same color or modify the color of a cluster for more appealing visualization. In addition, users can change parameter settings by dragging scroll bars provided in the interface in order to obtain promising result quickly.

4. Results

Our method was implemented on a PC with 3.50 GHz Intel Xeon E5-1620, 8 G memory, and an NVIDIA Quadro K4200. The GPU-accelerated ray-casting [30] volume algorithm was developed to render the results using C++ and GLSL shading language. Experiments were performed on various datasets to demonstrate the effectiveness of the proposed method. The datasets were obtained from the volume library (http://www9.informatik.uni-erlangen.de/External/vollib/).

4.1. Comparison

We first compared the proposed method with commonly used transfer function schemes, including basic intensity based method, traditional IGM based method using widgets, and a state-of-the-art IGM based method leveraging multilevel segmentation [7]. The experiments were performed on the Visible Human Male Head dataset. The results are shown in Figure 3. Figure 3(a) shows the visualization using basic intensity based method. It is observed that the bone and other interior organs are occluded by the skin. In addition, users have to extensively adjust TFs in a time-consuming trial-and-error way in order to achieve desired results. Figure 3(b) shows the visualization using traditional IGM based method, which demonstrates the advantage of using intensity-gradient histogram. For example, the left of Figure 3(b) is the initial result while the right of Figure 3(b) is the visualization by adjusting the sizes of the widgets and changing the opacity of each widget. This is a quite time-consuming task. Figure 3(c) shows the result using [7] which considers the intensity-gradient histogram as an image and segments it using graph cut to explore the different organs in the volume data. It is a hierarchical interaction process to remove the materials which users want to exclude in the visualization result. For example, as shown in Figure 3(c), if one wants to achieve the result in the right side, he (she) needs to perform at least six interactions on the initial result (left of Figure 3(c)). Figure 3(d) shows the visualization of our method, which can achieve satisfaction result automatically and users can further polish the result by performing very simple interactions.

4.2. Rendering Results
4.2.1. Visible Human Male Head

Figure 4 shows more visualization results of the Visible Human Male Head dataset () using our method. In order to more clearly show the effectiveness of our classification and interactivity approaches, we display each organ independently by setting the transparency of other classes to 0. Figure 4(a) shows the design of TFs and Figure 4(b) shows the visualization result of the whole volume. The red, green, yellow, light blue, and mazarine represent the skin, sinus, bones in the outside, bones in the inside, and teeth, while the grey represents the noise and nonsignificant tissues located at the boundaries of the volume. As teeth and the chin bone are quite close in the volume, our method classifies them into one cluster. Note that as the skins distance to volume data center is larger than other tissues and its gradient magnitude is also quite large, the opacity of the skin is set small in our method and thus it looks transparency in the visualization result to avoid occluding other inside organs.

4.2.2. Knee Dataset

Figure 5 shows the visualization of the knee dataset (). The dataset is acquired from CT scan and is used for anterior tibia osteotomy. The blue, red, and yellow in Figures 5(c), 5(d), and 5(e) represent the skin, tibia, and fibula, respectively, while the gray represents the structures excluded by the threshold.

4.2.3. Abdomen Dataset

Figure 6 shows the visualization of the abdomen dataset (). The dataset was acquired from CT scan of the abdomen and pelvis. A doctor must identify the key structures from it before performing surgical planning based on it. The blue, red, and green represent the skin, blood vessel, and pelvis, respectively. It is observed that the proposed method can automatically produce clear and meaningful visualization result for further diagnosis and treatment.

4.2.4. Other Datasets

In addition to CT scan data, our method is general enough to be employed for the visualization of other imaging modalities. Figure 7(a) shows the visualization of a computed tomography angiography (CTA) dataset (). Note that our method can visualize the small aneurysm in the middle of the brain (labeled by the green ellipse in the figure), which is quite important for diagnosis and treatment planning. Figure 7(b) shows the visualization of a MRI data.

The threshold of the spatial variance of a bin is a key parameter in the proposed method. It determines which bins should be eliminated from histogram. Figures 8 and 9 show the results of the Visible Human Male Head dataset and a foot dataset with different thresholds. The smaller the threshold is, the more the bins will be eliminated from the histogram. Figures 8(a)8(d) show the visualization results of the Visible Human Male Head dataset when we set the threshold as , , , and , respectively. Figures 9(a)9(c) show the visualization results of the foot dataset when we set the threshold as 0.75, 0.33, and 0.21, respectively.

4.3. Performance

Six datasets have been tested using our transfer function and its effectiveness in improving rendering quality has been demonstrated in Figures 39. To evaluate our method’s effectiveness of real time, we do the time performance test on these datasets. We also evaluate the time performance of the proposed method on these datasets. The results are shown in Table 1 with some key parameters. The time of computing in (2) is also tested as the preprocessing time.

5. Discussion

As we mentioned in Section 2, some previous works have been proposed based on IGM histogram to help transfer function design. Figure 3 shows the comparison of our method and most common approaches based on IGM histogram. From this comparison, we can see the traditional IGM based method needs to spend a lot of time to achieve satisfactory results while this is not convenient in clinical practice. Particularly in Figure 3(b), it is difficult for radiologists and physicians to find an appropriate location to put the widget, as they usually have no good understanding of the histogram. Figure 3(c) also needs a lot of interactions to remove the tissues which users may not want to observe. In our approach, it can remove noisy voxels by presetting the threshold of the spatial variance of a bin and the classification of voxels is automatically by using AP clustering that is shown in Figure 3(d). In Figure 6, the blood vessel is so small in intensity-gradient histogram that the users cannot easily put widget using traditional IGM based approaches, while our method can automatically figure it out using combined IGM and spatial information. This is convenient for doctors to operate to further improve the efficiency of diagnosis. Our method can be applied in many clinical applications such as maxillofacial surgery [31], tumor detection [32], and wearable device design [33].

Except the convenient operation, the other advantage of our method is that the spatial information of each voxel is also considered as the measurement of similarity and this can help distinguish different tissues, especially for the tissues whose intensity and gradient magnitude values are similar, such as bone and teeth that are shown in Figures 4(f) and 4(g). In Figure 5, the segmentation of tibia and fibula is a quite challenging task as they belong to the same material, while in our method by taking both intensity and gradient information and spatial information into consideration, the tibia and fibula can be successfully segmented and hence the doctors and acquire more information from the visualization result for diagnosis and treatment planning.

Although it needs some time to compute the spatial variance of a bin, our time performance test shows that the preprocessing time cost is so little that does not affect the real time rendering, and the variance of each bin is stored in a file when they computed firstly to avoid repeated computing.

6. Conclusion

In this paper, we present a novel automatic transfer function design scheme for medical data visualization based on the IGM histogram. Compared with previous studies also based on IGM histogram, our method considers both the intensity-gradient information and the spatial information of voxels, making the clustering of the voxels more accurate and effective. The AP clustering algorithm is employed to automatically generate the TFs. Compared with previous clustering algorithms employed in volume visualization, the AP clustering algorithm has much faster convergence speed and can achieve more accurate clustering results. To achieve meaningful clustering results, two novel similarity measurements, IGM similarity and spatial similarity, are proposed and integrated into the AP clustering algorithm. Experiments demonstrate the proposed method enable users to efficiently explore volumetric medical data with minimum interactions. Future investigations include assessing our method on more clinical data and attempting to automatically optimize the parameters based on the analysis of the input volume data.

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


The work described in this paper was supported by a grant from the Shenzhen-Hong Kong Innovation Circle Funding Program (nos. SGLH20131010151755080 and GHP/002/13SZ), a grant from the National Natural Science Foundation of China (Project nos. 61233012 and 61305097), a grant from the Guangdong Natural Science Foundation (no. 2016A030313047), and a grant from Shenzhen Basic Research Program (Project no. JCYJ20150525092940988).