Flow diagram
The flow chart of the proposed method is shown in Fig. 1. It includes the Gabor filter design, detection of spinal curves, localization of IVDs, and segmentation of IVDs. First, a set of 2D-Gabor filters, with different frequencies and directions, are used in image filtering to get a series of Gabor images. Second, Gabor features of the spine are calculated and the spinal curves are detected. Third, the Gabor features images (GFI) of IVDs are calculated and limited by spinal curves. Fourth, the IVD localization is performed by cluster analysis. Correction of localization is done to improve the accuracy of localization results. Last, segmentation of IVDs is performed based on the Gabor filter, localization results, and an improved self-adaptive threshold.
Setting of the Gabor filter
Effective extraction of the IVD features depends on the setting of the Gabor kernel function. Due to the elliptic shape of IVDs, the parameter \(\sigma\) of Gabor kernel is different in horizontal \(x\) and vertical \(y\) direction. In addition, there is an angle between the IVDs and the horizontal plane. Therefore, the Gabor kernel is defined as [25]:
$$\begin{aligned}\psi (x,y,\theta_{\mu } ,\omega_{\nu } ) &= \frac{1}{{2\pi \sigma_{x} \sigma_{y} }}\exp \left\{ { - \frac{1}{2}\left[ {\left(\frac{{x^{\prime}}}{{\sigma_{x} }}\right)^{2} + \left(\frac{{y^{\prime}}}{{\sigma_{y} }}\right)^{2} } \right] + i\omega_{\nu } x^{\prime}} \right\}\\ & \qquad \mu = 0, \ldots ,S - 1, \quad v = 0, \ldots ,K - 1 \end{aligned}$$
(1)
where, \(x^{\prime} = x\cos \theta_{\mu } + y\sin \theta_{\mu } ,y^{\prime} = - x\sin \theta_{\mu } + y\cos \theta_{\mu }\) represent the spatial locations of pixels. In spatial domain, the parameters \(\theta_{\mu }\), \(\omega_{\nu }\), and \(\sigma\) represent the direction, wavelength, and Gaussian window of Gabor filter bank, respectively.
The parameters of Gabor are set as follows: to describe the local features of images, a Gabor filter bank has \(S\) directions and \(K\) scales (usually \(S = 8,\,\,K = 5\)) [20]. However, the IVDs have a relatively uniform size and the differences among their angles are smaller. A Gabor filter bank in 16 directions with five scales is used. The upper frequency limit \(\omega_{\text {max} } = \frac{\pi }{2}\) is set according to the experience. Since the size of IVDs is relatively uniform, the frequency spacing factor \(f = \sqrt[4]{2}\) is set to obtain the effective center frequency. The effective radius of the Gaussian window is expressed as \(r_{v} = \frac{2\sqrt 2 \sigma }{{\omega_{v} }}\). According to the characteristics of the IVD (the width occupies at least 25 pixels and the thickness is about half of the width), \(\sigma_{x} = \frac{3k}{{\omega_{v} }},\sigma_{y} = \frac{6k}{{\omega_{v} }}\) (where \(k = \sqrt {2\ln 2}\)). The symmetric Gabor kernel window with the size of \(31 \times 31\) is used, based on experiment. Gabor filter bank is shown as Fig. 2.
Localization of IVDs
Detection of spinal curves
In order to reduce the impact of background on localization and segmentation, the spinal edges are detected to narrow the searching range. Since the spinal curves are nearly vertical, the spinal GFI is obtained by subtracting the GFI of horizontal direction from that of vertical direction. Taking \(G_{v,\mu }\) as the Gabor-filtered image with direction \(\mu\) and scale \(v\), the spinal GFI \(G_{spine}\) can be described as Eq. (2):
$$G_{spine} = \overline{{\sum\limits_{{\nu \in C,\mu \in U_{1} }} {G_{\nu ,\mu } } }} - \overline{{\sum\limits_{{\nu \in C,\mu \in U_{2} }} {G_{\nu ,\mu } } }}$$
(2)
where \(v \in C = \left\{ {0,1, \ldots ,S - 1} \right\}\), \(\mu \in U = \left\{ {0,1, \ldots ,K - 1} \right\}\), \(U_{1} = \left\{ {0,1,2,3,4,5,14,15} \right\}\), and \(U_{2} = \left\{ {7,8,9,10,11} \right\}\) represent the directions close to the vertical and horizontal planes, respectively. Obviously, the negative values in \(G_{spine}\) are not the edge information of spine. Then, the \(G_{spine}\) is processed as Eq. (3).
$$G_{spine} (x,y) = \left\{ \begin{array}{ll} G_{spine} ( x,y),&\quad G_{spine} ( x,y ) \ge 0 \\ 0,&\quad G_{spine} ( x,y ) < 0 \\ \end{array} \right.$$
(3)
$${G\left( n \right) = \sum\limits_{x = 1}^{n} {\sum\limits_{y = (M - p)/2}^{(M + p)/2} {G_{spine} (x,y)} } ,} \quad {n = 1 \ldots } N$$
(4)
where \(G(n)\) is the sum of the GFI for the first \(n\) column (Fig. 3b), and \(N\) is the column number of the image. Since the middle of the spine is nearly straight, the \(p\) rows in the middle of the MRI images are selected to detect the spinal edges. As shown in Fig. 3b, the curve of \(G(n)\) is approximately in the horizontal plane in the middle of the spinal region (the middle of the two white lines in Fig. 3a). The process include: first, the centers of the left and right spinal edges (the white crosses in Fig. 3a) are obtained by searching the crossing points for \(G_{spine}\) and center of spine (the white point in Fig. 3a from center to both sides). Second, the search range is limited between the two white lines in Fig. 3a. Third, the left and right spinal edges are obtained by searching the corresponding region from above the centers of left or right spinal edges to both sides, respectively. The Fig. 3c shows the result of spinal edges.
Localization of IVDs
The localization of IVDs is based on GFI information. GFI for localization of IVDs is calculated similar to the process for spines, but the directions of \(U_{1}\) and \(U_{2}\) are different. According to the anatomical knowledge, the angle between the long axis of the IVD and the horizontal line of MRI is almost ± 30°. Therefore, \(U_{1} = \left\{ {7,8,9,10,11} \right\}\) and \(U_{2} = \left\{ {1,2,3,4,5} \right\}\) represent directions close to the long axis of IVDs and close to the vertical directions, respectively. The IVDs information image \(G_{Mdisc}\) (Fig. 4b) is obtained by median filtering (elliptical filter template, long axis: 44 pixels, minor axis: 17 pixels) on the IVDs GFI (Fig. 4a). The searching range of IVDs is narrowed by spinal curves.
The process of localization needs some a priori information. The centroid of the IVD is the locating point. In order to obtain the coordinate value of the locating point \((X_{IVD} ,Y_{IVD} )\), two aspects of a priori anatomical knowledge are adopted: I)Abscissa \(x\) of most IVDs is roughly located in the same vertical line, except for the last few ones that lean slightly to the right; II) The offsets of vertical \(y\) range from 20 mm to 50 mm (25 pixels to 60 pixels in this study). The information of local maximum of \(G_{Mdisc}\), combined with a priori knowledge, is used for coarse localization. In this process, correction of centre-of-gravity shift is also performed to improve the localization accuracy in case of severe spinal curvature. The process includes:
-
Step 1:
Calculate the horizontal cumulative curve \(G_{h} \left( n \right)\) and the candidate vertical coordinate values. \(G_{h} \left( n \right)\) is obtained by adding each row with Eq. (5), Where \(M\) is the number of rows of the image. The coordinates of local maximum values are the candidates.
$${G_{h} \left( n \right) = \sum\limits_{x = 1}^{N} {G_{Mdisc} \left( {x,n} \right),} } \quad {n = 1 \ldots M}$$
(5)
-
Step 2:
Calculate the coarse \(Y_{IVD}\). According to the a priori information II, the closer points merge into one. The \(Y_{IVD}\) is achieved after removing the points which do not meet the a priori information II.
-
Step 3:
Calculate the coarse \(X_{IVD}\). The horizontal coordinate \(X_{IVD}\) is calculated by a process similar to that for \(Y_{IVD}\), but the regions of calculation are different. The mid-values of adjacent two \(Y_{IVD}\) are taken as boundaries to intercept the IVD part. The same action is performed as Eq. (6) in columns between the two boundaries. The \(X_{IVD}\) is obtained based on the a priori information I.
$${G_{v} \left( n \right) = \sum\limits_{y = 1}^{M} {G_{Mdisc} \left( {n,y} \right),} } \quad {n = 1 \ldots N}$$
(6)
-
Step 4:
Calculate the boxes of IVDs to correct the coarse(\(X_{IVD}\), \(Y_{IVD}\)). Based on the results of Step 3, the rectangular regions of IVDs (Fig. 4d) are found by the boundary of the local peak of \(G_{h} \left( n \right)\) (Fig. 4c) and \(G_{v} \left( n \right)\). The upwards and downwards edges of the box in Fig. 4c are calculated by the minimum points and zero points of \(G_{h} \left( n \right)\). Then the left and right edges of the box are calculated by the minimum points and zero points of \(G_{v} \left( n \right)\) limited by upwards and downwards edges. Figure 4c shows the boxes in \(G_{Mdisc}\). Figure 4d shows the boxes in the original image. The angles of IVDs are obtained by analyzing the GFIs information in the rectangular regions.
-
Step 5:
Calculate the accurate center of IVDs \((X_{IVD} ,Y_{IVD} )\). The binary image of \(G_{Mdisc}\) is calculated in the boxes obtained in step 3. Figure 4e shows the superposition of the original image and the binary image. Subsequently, the corrected localization of IVDs is obtained by calculating the centroid of the binary image. As shown in Fig. 4f, the white circles are the final localization results \((X_{IVD} ,Y_{IVD} )\).
Segmentation of IVDs
Combining the localization results and spinal curves,the segmentation of IVDs is performed based on GFIs with an adaptive threshold. The main steps include:
-
Step 1:
Delineate the candidate region of the IVD. The initial candidate region of each IVD for localization is comprised of the spine curves and the rectangular regions of IVDs (Fig. 4d). In this region, the average of GFIs of different frequencies \(\bar{G}_{\mu }\) is calculated in the direction \(\mu\). The maximum \(\bar{G}_{\mu }\) can reflect the boundary of the IVD. The \(\bar{G}_{\mu }\) in the direction \(\mu\) with \(S\) different frequencies is defined as Eq. (7).
$$\bar{G}_{\mu } \left( {x,y} \right) = \frac{1}{S}\sum\limits_{v \in C} {G_{\nu ,\mu } }$$
(7)
The binary image of maximum \(\bar{G}_{\mu }\) is obtained using the Otsu method. The up and down boundaries of the candidate region are comprised of the spine curves. The left and right boundaries are comprised of the fitting curve of the outer boundaries of the binary image of maximum \(\bar{G}_{\mu }\).
-
Step 2:
Calculate the adaptive local threshold \(T_{IVD}\). Due to the ambiguous boundary and diverse shapes of IVDs, global threshold is not appropriate to segment the IVDs. Therefore, the local threshold (\(T_{IVD}\)) of each IVD is adopted to segment the IVDs. As shown in Fig. 5, \(T_{IVD}\) is obtained by self-adaption iteration and used for segmentation of IVDs. The initial threshold is calculated by the Otsu method. \(A_{1}\) is the difference between the area of candidate region of the IVD and the area of the binary image of maximum \(\bar{G}_{\mu }\). \(A_{2}\) is 1/2 area of the binary image of maximum \(\bar{G}_{\mu }\). \(T_{1}\) and \(T_{2}\) are \(\hbox{min} (A_{1} ,A_{2} )\) and \(\hbox{max} (A_{1} ,A_{2} )\) area, respectively. The number of iterations is 30.
-
Step 3:
Calculate the coarse segmentation of the IVD. The coarse segmentation results of IVDs are used in steps 2 to 5, in Fig. 5, with an adaptive local threshold \(T_{IVD}\).
-
Step 4:
Post-process for IVD segmentation. The coarse segmentation result may contain holes, non-smooth areas, and superfluous portions. In order to obtain the accurate segmentation, a morphological operator is used: holes filling, erosion, dilation, and five foreground pixels on the background pixel indicate that the pixel is the foreground pixel. When a segmentation result has many connected regions, the largest region is retained.