3D Lumbar Intervertebral Disc Segmentation from MRI Data Sets

mm and the data set size is $40\times 512\times 512$ . The advantage of working with such datasets is that different channels provide complementary information for our disc segmentation task. In our proposed segmentation strategy, we always first extract either intensity or feature information about different tissues on each channel and then combine the 4 channel data into a single dataset.

Fig. 1

The four aligned channels of a patient data (for visualization purpose, we only show the middle sagittal (mid-sagittal) slice of each channel.)

Fig. 2

Initialization and 2D lumbar spine column detection. a User initialization by picking two landmarks indicating the centers of L1 and L5 in the middle sagittal slice. b Probability assignment (displayed as grey values) of the bone tissue in the mid-sagittal slice for 2D lumbar spine column detection. c 2D lumbar spine column detection result using the graphical model based detection algorithm, blue and green rectangles representing the vertebral bodies and IVDs respectively. d Graphical model of 2D lumbar column detection (Color figure online)

2.2 Lumbar Spine Column Identification

On the mid-sagittal slice, two landmarks are picked to indicate the centers of L1 and L5 vertebral bodies as shown in Fig. 2a. Starting from the initialization, we first carry out a 2D vertebral body and disc identification to localize vertebrae L1-L5 and the 5 target discs from the mid-sagittal slice. The geometrical information of the 2D identification is then used to guide a further 3D lumbar spine column modeling.

2.2.1 2D Vertebral Body and Disc Identification

Solutions for spine location and disc labeling include feature-based bottom-up methods, statistical model-based methods and graphical model-based solutions. For a detailed review of the existing methods, we refer to [14]. In this paper, the 2D vertebral body and disc identification is achieved using a graphical model based strategy introduced in [14]. Compared with the graphical models in [5, 6], the advantage of the graphical model in [14] is that both the low level image observation model and the high level vertebra context potentials need not to be learned from training data. Instead they are capable of self-learning from the image data during the inference procedure. For completeness, here we describe the key components of the method that we previously introduced [14].

The graphical model: The graphical model is given in Fig. 2d. Each node $V_{i}$ represents a connected disc-vertebrae-disc chain of the spine, whose geometrical parameters are given by $X_{i}$ . We define

The component observation model. $p(I|X_{i})$ of a single component $V_{i}$ representing the probability that the configuration $X_{i}$ of the node $V_{i}$ match the observed images .
The potentials. $p(X_{i},X_{j})$ between neighboring components $V_{i}$ and $V_{j}$ encoding the geometrical constraints between components which are defined by the anatomical structure of the spine column.

The identification of the spine column from the mid-sagittal slice can then be formalized as to find the optimal configurations of $\{V_{i}\}$ , $X=\{X_{i}\}$ that maximize

$\begin{aligned} P(X|I)\propto \varPi _{i}p(I|X_{i})\varPi _{i,j}p(X_{i},X_{j}) \end{aligned}$

(1)

with

$\begin{aligned} p(I|X_{i})=p_{I}(I|X_{i})p_{G}(I|X_{i}) \end{aligned}$

(2)

and

$\begin{aligned} p(X_{i},X_{j})=p_{S}(X_{i},X_{j})p_{O}(X_{i},X_{j})p_{D}(X_{i},X_{j}) \end{aligned}$

(3)

$p_{I}(I|X_{i})$ and $p_{G}(I|X_{i})$ stand for the probabilities that the observed image intensity and image gradient distributions match the geometrical parameters $X_{i}$ respectively. $p_{S}(X_{i},X_{j})$ , $p_{O}(X_{i},X_{j})$ and $p_{D}(X_{i},X_{j})$ are the geometrical constraints on the sizes, orientations and distances between neighboring components. All the observation models and constraints can be designed according to the observed data and prior anatomical knowledge of the spine structure. For detailed formulation of these terms, we refer to [14].

Optimization: The optimization is achieved as an inference on the graphical model. The method introduced in [14], which is a particle based nonparametric belief propagation on the graphical model, is used here to carry out the inference.

Figure 2b shows an example of the bone tissue probability assignment on the mid-sagittal slice, which is computed from the user supplied 2 landmarks (Fig. 2a) and 4 channel volume data (see Fig. 1 for an example) by a Gaussian distribution modeling and an equally weighted combination of the intensity distributions of the bony tissue in the 4 channels. This image is used for the computation of the intensity observation model $p_{I}(I|X_{i})$ during the 2D lumbar column detection. Figure 2c gives the 2D lumbar column detection result. It can be observed that the centers, sizes and orientations of the vertebral bodies and IVDs are correctly identified.

2.2.2 3D Lumbar Spine Column Modeling

We model each lumbar vertebral body as an elliptical cylinder and the lumbar spine column as a variable-radius soft tube. Details of the modeling procedure are described as follows:

3D modeling of each vertebral body: From the 2D vertebral body identification results, the position, hight, radius and orientation of each vertebral body and the image intensity distribution of the bone region can be estimated by modeling the vertebral body as a cylinder. Accordingly for each voxel in the neighbourhood of the estimated cylinder, we can assign the probability that it belongs to the bony tissue. To refine the 3D modeling of the vertebral bodies, we then further model the vertebral body as an elliptical cylinder, a least-squares geometric fitting to the voxels assigned with a high probability (3a.

Given the 3D soft-tube lumbar spine column model, the spine column region can be extracted from the observed data sets (Fig. 3b). By further eliminating the bony tissue region using the 3D models of vertebral bodies, the candidate region for each target disc can be localized as shown in Fig. 3c–f. The following 3D IVD segmentation is then carried out on the extracted candidate IVD regions.

Fig. 3

3D lumbar spine column detection and modeling. a the 3D soft-tube model of the lumbar spine column; b segmented lumbar spine column image; c–d segmented disc candidate regions in sagittal slices; e–f segmented disc candidate regions in coronal slices. Although all tasks are conducted in 3D, here we show the results in 2D slices for visualization purpose