Robust MRI abnormality detection using background noise removal with polyfit surface evolution

Image segmentation plays a vital role in MRI abnormality detection. This paper presents a robust MRI segmentation method to outline potential abnormality blobs. Thresholding and boundary tracing strategies are employed to remove background noises, and hence, the ROIs in the whole process are set. Subsequently, a polyfit surface evolution is proposed to approximately estimate bias field, which makes segmentation robust to image noises. Simultaneously, customized initial level set functions are devised so as to detect subtle bright and dark blobs which are highly potential abnormality regions. The proposed method improves bias field estimation and level set method to acquire fine segmentation with low computational complexities. Analysis of experimental results and comparisons with existing algorithms demonstrates that the proposed method can segment weak-edged, low-resolution MR brain images, and its performance prevails in accuracy and effectiveness.

With the exploration of magnetic resonance imaging (MRI) technology, dramatic changes have been made in visualization of anatomical structures [1, 2]. MRI is a powerful imaging modality driven by its flexibility and sensitivity to a broad range of tissue properties [3]. However, a CT scan is best suited for viewing bone injuries, diagnosing lung and chest problems, and detecting cancers [4, 5]. MRI is used to find other problems, such as tumors, bleeding, injury, blood vessel diseases, and infections. For example, researchers strive to detect brain abnormalities in MR images [6, 7]. In this paper, our goal is to segment MR brain images to assist doctors’ diagnosis.

Intensity inhomogeneity or bias field in MR image, which arises from the imperfections of the image acquisition process and manifests itself as slow intensity variation over the image domain [8]. This inherent artifact is difficult for a human to perceive. However, many segmentation methods are very sensitive to spurious intensity variations. In 1986, Haselgrove and Prammer [9] first discussed MRI intensity inhomogeneity correction. However, Haselgrove and Prammer assumed that the position and orientation of the surface coil were known, as obstructed its practical application. In the past two decades, many bias field correction algorithms have been presented. The authors in [8] noted that two main approaches, namely prospective [10–14] and retrospective approaches [15–21], had been applied to minimize the intensity inhomogeneity in MR images. Among these abovementioned methods, Andersen et al. [15] incorporated radio frequency inhomogeneity correction into probabilistic classification model to partition an MR image. In addition, Li et al. [22] proposed a new approach for bias field estimation and tissue segmentation in an energy minimization framework. They mathematically justified that the proposed energy is convex in each of its variables. Quantitative evaluations and comparisons verified its robustness and accuracy. The bias field estimation proposed by Li et al. needs no previous knowledge of noise distribution, which makes it more applicable. However, it is a region-based segmentation method and cannot extract small homogeneous regions as level set methods do.

Image segmentation is a fundamental process in many medical imaging applications [23]. Segmentation methods can be roughly divided into eight categories [24–26]: (a) thesholding approaches, (b) region growing approaches, (c) classifiers, (d) clustering approaches, (e) Markov random field models, (f) artificial neural networks, (g) deformable models, and (h) atlas-guided approaches. The level set method, as one of the deformable models, was first introduced to delineate region boundaries by using closed parametric curves [27]. It is able to acquire closed contours of regions from an image, which helps partition of a medical image accurately. The authors in [28–30] considered a two-phase level set formulation, in which only one level set function was used to construct two membership functions to segment the image domain into two disjoint regions. The two-phase level set method can only partition images into two parts, making it unsuitable for multi-class segmentation in some medical applications. Recently, researchers have developed multi-phase level set methods [31, 32], using two or more level set functions to define more than two membership functions, which makes it more practical in medical applications. Li et al. [31] have combined bias field estimation with multi-phase level set method. In this framework proposed by Li’s, the formation of initial level set functions was not mentioned yet. Experimental results show initial level set functions influence final segmentation to some extent, especially in certain small regions. At the same time, bias field estimation in [31] lacks denoising capability as well.

Focusing on noise depression and widely applicable initial level set functions, we import ROI (region of interest), bias field polyfit estimation, and customize initial level set functions to multi-phase level set methods. The main contributions in this paper include:

1. (1)

   A coarse to fine segmentation method is proposed. Otsu’s method is employed not only to remove background noise but also to set the ROI for the level set method.

2. (2)

   A multi-phase level set method with bias field polyfit estimation is proposed. This method can partly depress image noise latent in tissues.

3. (3)

   Customized initial level set functions are formulated. The initial level set functions are widely applicable to detect small tissues in MR images.

This remainder of this paper is organized as follows. We first review MRI bias field estimation and multi-phase level set method in Section 2. Section 3 presents a robust MRI abnormality detection method, including background removal or ROI setting, multi-phase level set method with bias field polyfit estimation, and customized initial level set functions. Experimental results and comparison with existing methods are outlined in Section 4. Finally, concluding remarks are given in Section 5.

In previous works [8, 22, 31], the model of an MR image has been widely accepted as:

where b(x) is bias field that accounts for the intensity inhomogeneity in the observed image I(x), n(x) is the additive Gaussian noise with zero mean, J(x) is the true image without bias field and noise, and x is a tuple like (x,y) representing image coordinates.

An energy minimization formulation [22] in the image domain Ω for bias field estimation is defined by:

In this paper, the domain of integration Ω is omitted for simplicity. Obviously, the solution for the multiplicative intrinsic components b and J is an ill-posed problem without constraints on the variables b and J. Thus, Li et al. [22] have made the two following assumptions:

1. (1)

   Image domain Ω can be partitioned into disjoint parts as described below:

   $$\Omega = \sum \limits_{i=1}^{N} \Omega_{i},~\Omega_{i} \bigcap \Omega_{j} = \Phi~(i \neq j), $$

   in which J is approximately constant c _i for the ith region Ω _i . Furthermore, c _i is regarded as class center of fuzzy c-means algorithm [33].

2. (2)

   The bias field b varies slightly in a small neighborhood region of any point y on the ith region Ω _i , and leads to the approximate relation in Eq. (2):

   $$b(\mathbf{x}) J(\mathbf{x}) \approx b(\mathbf{y}) c_{i},~\mathbf{x} \in \mathcal{O}_{\mathbf{y}} \bigcap \Omega_{i} $$

   where \(\mathcal {O}_{\mathbf {y}} \triangleq \{ \mathbf {x}:|\mathbf {x} - \mathbf {y}| \leq \rho \} \).

Furthermore, a non-negative window function K(y−x) is introduced, being K(y−x)=0 for \(\mathbf {x} \notin \mathcal {O}_{\mathbf {y}}\), to represent the energy function bounded on region Ω _i as follows:

Finally, extending the integral domain y∈Ω _i in (3) to the entire image domain, the energy minimization in (2) is formulated as:

Li et al. [31] proposed level set functions to construct membership functions, which form partitions Ω _i . For the case of N ≥ 3, two or more level set functions ϕ ₁,…,ϕ _k are used to define membership functions M _i , given by:

Putting (5) into (4), the energy functional is given by:

where Φ=(ϕ ₁,…,ϕ _k ) is a vector valued function.

Reversing the order of integration in (6), the energy functional can be rewritten as:

Denoting the constants c ₁,…,c _N by a vector c=(c ₁,…,c _N ), Eq. (7) is simplified as:

where \(e_{i}(\mathbf {x}) = \int K(\mathbf {y}-\mathbf {x}) |I(\mathbf {x}) - b(\mathbf {y}) c_{i}|^{2} d\mathbf {y}\).

The level set formulation is given as follows:

The energy term L(ϕ _j ) in (9) computes the length of the zero level set of ϕ _j in the conformal metric and it is defined by:

where ∇ is the vector differential operator, namely gradient operator, H is the heaviside function.

The energy term R(ϕ _j ) in (9) is introduced in [28] to formulate level set evolution without re-initializations. R(ϕ _j ) is given by:

The energy functional \(\mathcal {F}(\boldsymbol {\Phi },\mathbf {c},b)\) is minimized iteratively with respect to each of the variables Φ,c,b, given that the other two have been updated in the previous iteration.

1. (1)

   Optimization of ϕ _j :

   $$ \left\{ \begin{array}{l} \frac{\partial \phi_{j}}{\partial t} = - \sum \limits_{i = 1}^{N} \frac{\partial M_{i}(\boldsymbol{\Phi})}{\partial \phi_{j}} e_{i} + \nu \delta(\phi_{j}) \text{div} \left(\frac{\nabla \phi_{j}}{| \nabla \phi_{j} |} \right)\\ + \mu \left[ \Delta \phi_{j} - \text{div} \left(\frac{\nabla \phi_{j}}{|\nabla \phi_{j}|}\right) \right]\\ e_{i}(\mathbf{x}) = I^{2} \cdot (\mathbf{1} \ast K) - 2c_{i} I \cdot (b \ast K) + c_{i}^{2} (b^{2} \ast K) \\ (j = 1,2,\ldots,k) \end{array} \right. $$

   (12)

   where · calculates the per-element product of two matrices, ∗ is the convolution operator, 1 is a matrix of ones with the same size as I, δ(·) is the derivative of the heaviside function and div(·) is the divergence.

2. (2)

   Optimization of c:

   $$ c_{i} = \frac{\int (b \ast K) I M_{i}(\boldsymbol{\Phi} (\mathbf{y})) d \mathbf{y}}{\int (b^{2} \ast K) M_{i}(\boldsymbol{\Phi}(\mathbf{y})) d \mathbf{y}} ~(i = 1,\ldots,N) $$

   (13)
3. (3)

   Optimization of b:

   $$ b = \frac{IJ^{(1)} \ast K}{J^{(2)} \ast K} $$

   (14)

   where

   $$\left\{ \begin{array}{ll} J^{(1)} = \sum \limits_{i=1}^{N}c_{i} M_{i}(\boldsymbol{\Phi}(\mathbf{y}))\\ J^{(2)} = \sum \limits_{i=1}^{N}{c_{i}}^{2} M_{i}(\boldsymbol{\Phi}(\mathbf{y})) \end{array} \right.. $$

Repeated trials expose imperfections of Li et al.’s method [22, 31] in detecting abnormalities from MR images: (a) background noise gives rise to redundant contours, (b) unsmooth bias field estimation exists, and (c) some possible abnormalities are omitted. We propose a robust MRI abnormality detection method to fix these problems based on a multi-phase level set method in combination with bias field estimation.

3.1 Background removal (ROI setting)

As we know, MRI is vulnerable to noise contamination. The MRI signal is complex valued, with the real and imaginary components affected by independent Gaussian noises. The magnitude of the noise of this complex valued signal has a Rician distribution [34]. Assuming the signal in the image background is zero, Santiago et al. [35] estimated noise in a MRI using background pixels. If noise in background pixels is severe, it will influence final curve evolution using Li et al.’s method mentioned in Section 2, see Fig. 1 a.

Inaccurate curve evolution by Li et al.’s [31] method with background noises: a Final contours without background noise suppression. b Contours after background noise removal. c Contours by our method

Full size image

Fortunately, the gray values of MRIs differ greatly in background and organs. Thus, thresholding is a simple but efficient method to remove the background from MRI. In our method, we first extract coarse regions of a target using Otsu’s method [36]. Subsequently, we retrieve connected components from the binary image and select the contour outlines in the image whose areas are greater than a predefined A value (step (iii) below). Finally, we fill the regions bounded by the preceding selected contours, which serve as the ROI in this paper. Assuming gray values of background are much lower than the foreground ones, the procedure can be outlined as follows:

1. 1.

   Compute the optimal threshold value using Otsu’s algorithm and denote it by T.

2. 2.

   Apply a fixed-level threshold to each array element of image f and obtain the binary image g:

   $$ g(i,j) = \left\{ \begin{array}{rl} 255, & f(i,j) \ge \eta T\\ 0, & \text{otherwise} \end{array} \right. $$

   (15)

   where constant 0<η≤1 can guarantee that the regions of foreground with lower gray values are not ignored. Based on empirical observations, this constant commonly takes values in the interval [0.7,0.9].

3. 3.

   Trace only the outer contours of image g. For the ith contour, if the area is greater than A, fill the area bounded by the contour with the gray value 255 and denote the corresponding mask image by B _i . The constant A can skip small regions incurred by noise. Actually, A can be easily defined by professionals’ prior knowledge about the area of the target region. The union set of such B _i , denoted by \(B~=~\bigcup B_{i}\), is the final ROI image, with the notation R in which the value 1 represents the points to be further processed:

   $$ R(i,j) = \left\{ \begin{array}{rl} 1, & B(i,j) = 255\\ 0, & \text{otherwise} \end{array} \right. $$

   (16)

Background noise is likely to bring out spurious contours by Li et al.’s method, shown in Fig. 1 a. After the R is calculated, background noise is removed with the values of background pixels assigned to zero. However, it is apparent that there are still some faulty contours with open contours, see Fig. 1 b. For this reason, it is wise to neglect the background region which causes abnormalities. Furthermore, only the pixels in these locations whose corresponding values are 1 in the ROI R are considered in our proposed method, not only to obtain accurate final contours (see Fig. 1 c) but also to speed up calculation.

3.2 Multi-phase level set method with bias field polyfit estimation

In Section 3.1, ROI is introduced to eliminate the negative effect of background noises. Noise in the target region of a MRI does harm to the estimation of bias field. As described in (7), local neighboring pixels are contributive with the definition of the window function K. Confined to a local domain, the b can be approximately represented as a linear combination of simple functions [22]:

Naturally, we introduce bias field polyfit estimation to the multi-phase level set method. Experiments in Section 4 on detecting abnormalities from MR images demonstrate that this can significantly improve segmentation results.

In mathematics, polynomial bases have been widely used. In this paper, we employ the basis of polynomials with degree p that fits the bias field. It is presented as:

The relationship between the degrees p and terms L is described in Table 1.

Substituting (17) into (8), the energy functional \(\mathcal {E}\) is converted into:

where \(\mathcal {E}_{\mathbf {x}} = \sum \limits _{i=1}^{N} \int K(\mathbf {y} - \mathbf {x}) |I(\mathbf {x}) - \mathbf {w}^{\mathrm {T}} G(\mathbf {y}) |^{2} M_{i}(\boldsymbol {\Phi }(x)) d \mathbf {y}\).

Therefore, \(\frac {\partial \mathcal {F}}{\partial w} = - 2 \mathbf {v} + 2 \mathbf {A} \mathbf {w}\), where

Applying the convolution operation, the equations above will be given by:

Note that \(\mathbf {v} \in \ \mathbb {R}^{L \times 1}\), \(\mathbf {A} \in \mathbb {R}^{L \times L}\). Let \(\frac {\partial \mathcal {F}}{\partial \mathbf {w}} = \mathbf {0}\), we get:

and,

Bias field polyfit estimation can retain the smoothness, as shown in Fig. 2. Compared to existing literature [22] and [31], one of the differences in our approach is the additional convolution operations in Eq. (21), which help localize and remove noise.

Bias field estimation by a Li’s method [22]; b the proposed method in this paper

Full size image

3.3 Customized initial level set functions to detect abnormalities

One of the disadvantages of traditional level set methods is that it requires manual placement of a closed curve near the desired boundary. Recently, researchers have employed random initial level set functions. Normally, once level set functions are stable, the process of evolution stops. Nevertheless, it is a challenge to jump out of a local extremum in (9). Moreover, some abnormalities occur as negligible patches, therefore no suitable initial level set functions can result in such abnormalities being detected. As a result, we introduce level set functions which intersect the target region in MRIs as much as possible. When level set functions during adjacent iterations are invariant, corresponding contours are the exact boundaries of target regions including abnormalities. These initial level set functions help us find out the abnormalities automatically, especially the negligible patches.

In order to segment more regions of interest, such as bright and dark blobs, we formulate initial level set functions so as to make initial contours scatter linearly, equally spaced in image domain. Classification of anatomical structure appearances in a MRI slice of the brain [37], suggests that the ROI mentioned in this paper can be classified into light, gray, and dark parts; namely, two level set functions constructing three membership functions (k=2,N=3). The initial contours are illustrated in Fig. 3. Here, red and blue contours in Fig. 3 b are a visual presentation of the level set functions \(\phi _{1}^{0}\) and \(\phi _{2}^{0}\) respectively, which guarantees that initial contours intersect the tissues extremely, merge or split to evolve to the final desirable contour.

Initial contour: a ROI-based initial contour and b partial enlarged view

Full size image

The pseudo-code for definition of \(\phi _{1}^{0}(x,y), \phi _{2}^{0}(x,y)\) is shown in Table 2, where the image dimension is W×H, ρ>0 is a constant with any value, and the ROI is bounded by a rectangle: {(x,y)|x _min≤x≤x _max,y _min≤y≤y _max}. In this case, the relationship between membership functions and level set functions is given by:

It is worth pointing out that the member functions mentioned meet 0≤M _i (x)≤1 and \(\sum \limits _{i = 1}^{3} M_{i}(\mathbf {x}) = 1\).

3.4 ROI-based numerical implementation

In order to remove the influence of noise in the background and reduce computational cost, we propose an ROI-based implementation for the aforementioned algorithm.

In numerical practice, the Dirac function δ(x) in (12) is smoothed as:

Iteration for ϕ _j is proposed by:

In this paper, the window function K is taken as a Gaussian function with a standard deviation σ, denoted by K _σ . Our proposed ROI-based numerical implementation is presented as follows:

1. 1.

   ROI mask. I=I·R.

2. 2.

   Initialization. Set the following parameters: N, η, p, σ, ν, μ, ε and τ. Initialize b and ϕ ₁,…,ϕ _k . Due to their invariability during iterative procedures, we prepare G, G G ^T, 1∗K _σ , (1∗K _σ )·I ², (G∗K _σ )·I and G G ^T∗K _σ .

3. 3.

   Iterative procedure. For the ith iteration, first update c via (13). Second, update Φ via (25), (12), and (26). Finally, update b via (21), (22), and (23).

4. 4.

   ROI mask. ϕ _j =ϕ _j ·R for next iteration.

5. 5.

   Redo step (iii) until i exceeds the predefined maximal iteration.

We analyzed 112 frames of MR brain images; detailed information can be seen in Table 3. The algorithm was implemented using the Microsoft Visual Studio development platform with the Open CV (Computer Vision) library, and Matlab. In this paper, we have used the same parameters listed in Table 4, except the specific statement. In the experiments, the number for the maximal iteration needed is 50. The energy objective function in (9) keeps descending, shown in Fig. 4. This shows that our proposed method has good convergence.

Energy functional values in the iterations

Full size image

4.1 Segmentation results

An overview of the proposed method for MRI segmentation is shown in Fig. 5. Here, the original image is shown in Fig. 5 a, whose coarse segmentation result, namely ROI, is shown in Fig. 5 b. After 10 iterations, the curve evolves as shown in Fig. 5 c. After 50 iterations, the final contour and segmentation results are illustrated in Fig. 5 d–f. Three split parts of segmentation results in Fig. 5 g–i demonstrate the capability of our method to segment the bright and dark blobs (high likelihood of potential abnormality) simultaneously. Here, in order to retain the running results by our scripts, contours in red and blue imply that the boundaries are induced by different level set functions, as mentioned in Section 3.3. Without loss of generality, we henceforth adopt red contours later in this paper.

Segmentation processes for the 88th frame: a the original MR image, b ROI, c contour after 10 iterations, d final contour after 50 iterations, f segmentation results, g one part segmented, h bright blob in MR image, and i dark blob in MR image

Full size image

Experimental results for the other frames are shown in Fig. 6. The relevant parameters are employed as shown in Table 4. The segmentation performance verifies that the listed parameters are image independent.

Segmentation processes for the other frames: a the 14th MR image, b–d segmented parts, e the 50th MR image, f–h segmented parts, i the 67th MR image, and j–l segmented parts

Full size image

Outside the data set mentioned above, we experiment on the MR brain images [38], which includes 60 images (width 256, height 256, bit depth 16). Taking the 6th image (see Fig. 7 a) as an example, our proposed method can segment the images into parts: skull (see Fig. 7 b), white matter (see Fig. 7 c, d), and other parts.

Segmentation processes for the other frames: a the 14th MR image, b–d segmented parts, e the 50th MR image, f–h segmented parts, i the 67th MR image, and j–l segmented parts

Full size image

4.2 Segmentation assessment

For objective evaluation of our method, we consider manual contours drawn by specialists as the ground truth. We introduce two metrics to segmentation assessment. One metric m ₁ is defined to estimate the difference between the ground truth and the contours extracted by our method. For one specific ground truth contour, we first select its nearest contour edge in contours from our proposed method, then we calculate the distance between the contour pair. The other metric m ₂ is to account for how many accurate possible abnormal blobs are detected.

Assuming ground truth contours are denoted by \(C_{g}~=~\left \{ C_{g}^{1}, C_{g}^{2}, \ldots, C_{g}^{n_{g}} \right \}\), and contours from the proposed method are denoted by \(C_{p} = \left \{ C_{p}^{1}, C_{p}^{2}, \ldots, C_{p}^{n_{p}} \right \}\).

For the ith contour, the corresponding best matching contour in the proposed method is defined by:

where \(d\left (p_{j},C_{p}^{k}\right)\) is the distance from point p _j to contour edge \(C_{p}^{k}\).

Therefore, the measurement depicting the difference of the contours representing the ground truth and the proposed method is given by:

where \(\left |C_{p}^{k_{i}}\right |\) is the number of elements in \(C_{p}^{k_{i}}\).

Assuming N _p is the number of potential abnormality blobs detected by the proposed method and N _g is the number of possible abnormality blobs provided by specialists, the metric m ₂ is defined by the ratio of N _p to N _g :

The smaller m ₁ or the bigger m ₂ is, the better the segmentation result is. Comparison results between Li’s method and our method using these metrics are shown in Table 5 and Fig. 8.

Comparison of errors in different contours

Full size image

4.3 Comparisons with existing methods

Without less noise affected MRI [39], our and Li’s method [31] give similar contours, shown in Fig. 9, which demonstrates the validity of our method. Our method prevails in two aspects: (1) the outer contour is more smooth than Li’s method, and (2) it detects more subtle bright or dark blobs which are potential abnormalities.

Comparison with the work [31]: a original MR image, b result via [31], and c result based on our method

Full size image

Compared to bias field estimation [22]-based segmentation (see Fig. 10 b), our method produces more accurate segmentation result with fine localization of regions of interest, as shown in Fig. 10 d. However, result based on [22] cannot give the contour directly and it misses the black blob. Furthermore, it cannot output accurate closed regions, such as the gray matter region (polygon bounded by a closed yellow curve in Fig. 10 d). In comparison with our method, traditional level set method [28] only extracts the outermost contours, as shown in Fig 10 c. Our method detects complete contours (see Fig. 10 d), in which the contour is highlighted in yellow, including missed detections in [28].

Comparison with the works [22, 28]: a the 87th MR image, b result via [22], c result via [28], and d result based on our method

Full size image

The authors in [31] take intensity inhomogeneities and level set method into account for image segmentation. The difference between ours and the one in [31] lies in that we minimize the w, while the latter minimizes the b directly. Therefore, our method can segment more regions of interest, including bright and dark blobs simultaneously, as shown in Fig. 11, which helps to improve diagnosis in medical applications.

Comparison with [31]: a the 84th MR image, b the result via [31], and c the result based on our method

Full size image

4.4 Discussion

In this paper, we introduced Otsu algorithm to remove background noise from original image. It is simple but effective in coarse segmentation. Such course segmentation helps to improve accuracy in final segmentation. Additionally, bias field polyfit estimation gives a more continuous bias field map, see Fig. 2 b. Furthermore, customized initial level set functions proposed in this paper helps to segment weak-edged, low-resolution MR brain images. In a sense, we incorporate bias field estimation into multi-phase level set method to acquire fine segmentation with low computational complexity. The proposed method can outline possible abnormalities. However, lesion recognition is to be determined in future work.

In this paper, we proposed a robust MRI abnormality detection method, which utilized background noise removal and polyfit level set method. Note that we improve the work in [22, 31], with our method being able to segment subtle bright or dark blobs automatically. Background removal helped us to exclude background noise, which also provided ROI to speed up successive processing. Polyfit level set method with customized initial level set functions in this paper does segment MRIs into individual parts, including tiny tissues. Experimental results demonstrate that our method has good performance in extracting all the blobs in MRIs, in which, there are potential lesions. But, as for each divided part, whether it is abnormal is yet to be determined. In future work, we will apply our technique on brain injury detection, such as for white matter injury detection.

The authors would like to thank NSERC, Canada, for their financial support of this research.

Funding

The work was supported in part by the Open Project of the Artificial Intelligence Key Laboratory of Sichuan Province under grant nos. 2014RZY02 and 2017RZJ03, the Open Project of Sichuan Province University Key Laboratory of Bridge Non-destruction Detecting and Engineering Computing under grant no. 2014QZY01, Natural Science Foundation of Sichuan University of Science and Engineering (SUSE) under grant nos. 2015RC08 and 2017RCL23, and Educational reform project of SUSE under grant no. JG-1707. The work was also supported in part by the Training Programs of Innovation and Entrepreneurship of Sichuan Province for Undergraduates under grant no. 201710622088. The work was also supported in part by Guangxi Key Laboratory of Cryptography and Information Security under grant no. GCIS201607.

Availability of data and materials

Not applicable.

Authors and Affiliations

1. Artificial Intelligence Key Laboratory of Sichuan Province, Sichuan Province University Key Laboratory of Bridge Non-destruction Detecting and Engineering Computing, Zigong, China

   Changjiang Liu & Jun Ye

2. Sichuan University of Science and Engineering, Zigong, 643000, China

   Changjiang Liu & Jun Ye

3. Department of Computing Science, University of Alberta, Edmonton, T6G 2H1, Canada

   Changjiang Liu, Irene Cheng & Anup Basu

4. Guangxi Key Laboratory of Cryptography and Information Security, Guilin, 541004, China

   Jun Ye

Contributions

CL and AB designed the algorithm, and CL wrote the paper and also performed the experiments. IC and AB reviewed and edited the manuscript. JY gave a critical suggestion on experimental section. All authors read and approved the manuscript.

Corresponding author

Correspondence to Changjiang Liu.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions