Vertebral Segmentation Using Part-Based Decomposition and Conditional Shape Models

Fig. 1

Examples of segmentations (in blue) obtained with local PDMs, showing suboptimal fitting in areas of complex geometry and high curvature on the spine (Color figure online)

Amongst existing techniques for vertebral image segmentation, statistical models of 3D shape [4] have been extensively used [2, 3, 5, 7, 8] due to their ability to build a shape prior from a representative training population. However, these methods consider at best a whole vertebra as the smallest unit for the construction of the point distribution models (PDMs). Due to the large variability of the vertebrae in particular at the processes and the generally small number of samples available for training, the obtained models are too constraining and not flexible enough to localize the fine details at areas of high curvatures.

In this paper, we present a new method for detailed modeling and segmentation of the vertebrae. The fundamental idea behind the proposed technique is to decompose each vertebra into a set of subparts based on their geometrical properties. By using a Vononoi [9] decomposition approach, we ensure that each of the main subparts of the vertebrae is well identified. With this approach, the shape constraints are effectively relaxed, allowing for an improved encoding of the fine details and shape variability at all the regions of the structures. Subsequently, in order to maintain the statistical coherence of the ensemble, conditional models are used to model the statistical inter-relationships between the different subparts. For shape reconstruction and segmentation, a robust model fitting procedure is introduced to exclude outlying inter-part relationships in the estimation of the shape parameters. The proposed technique is validated with a total of 30 spinal CT scans.

2 Method

The proposed framework consists of three main stages. Firstly, in Sect. 2.1, a subdivision of each vertebra into a number of subparts is proposed based on an approximate Voronoi region decomposition. Subsequently, the conditional models describing the statistical inter-relationships between the subparts are presented in Sect. 2.2. Finally, a model fitting approach based on all pair wise conditional models is introduced in Sect. 2.3, with the aim to estimate the shape parameters for each subpart robustly.

2.1 Vertebral Decomposition

Let us denote $\mathbf{x}=(\mathbf{x}_1,\dots ,\mathbf{x}_n)^T$ the landmark-based shape representation of each vertebra, where $$n$$

is the number of landmarks used to discretize the 3D shape. The aim of this section is to obtain a subdivision of $\mathbf{x}$ into $$K$$

sub-components.

We do this in this paper by using a polygon clustering algorithm described in [9], which provides a compact subdivision of the shape based on the concept of a Voronoi diagram [1]. Furthermore, let us denote $V=(C_1,\dots ,C_m)^T$ the triangulation of the shape $\mathbf{x}$ , where $$C_i$$

, $i=1, \dots , m$ , represents each face on the mesh, and $$E_j$$

the set of edges between all adjacent triangles. Given the centroids $$c_i$$

corresponding to each of the triangular triangles $$C_i$$

, the algorithm computes an approximation of a centroidal Voronoi diagram (CVD). The energy term to minimize is:

$\begin{aligned} F = \sum _{k = 1}^K \left( \sum _{i \in R_k} w_i||c_i - c_{R_k}||^2 \right) , \end{aligned}$

(1)

where

is a subset of $$V$$

(i.e. subpart of the shape), $c_{R_k}$ is the center of the region, and $$w_i$$

is the area of triangle $$C_i$$

To minimize Eq. 1, we use an iterative approach over the subset of edges $$E_j$$

between adjacent regions. The regions $$R_k$$

are initialized as a single triangle that is randomly chosen amongst $$V$$

. The remaining triangles are assigned to the $$null$$

region

. We then iterate only over the edges that are between two regions $$R_k$$

and

, or between any $$R_k$$

and

. Then we assign one of the two triangles adjacent to the current edge to the region that minimizes Eq. 1. At some point, the region $$R_0$$

will become empty. The iterative algorithm will continue until no modification of the region assignments for the triangles lead to an improvement of the subdivision according to Eq. 1.

Given that the regions of highest complexity and curvature on the vertebra are the different vertebral processes, we chose the lowest number of clusters $$K$$

such that points selected roughly at the distal region of every process (transverse, superior/inferior articular and spinous) are assigned to a different region after the subdivision.

Once the algorithm converges, all $C_i \in V$ will belong to a unique region $$R_i$$

. However, the points that lie on the boundary edges between subparts will now belong to two adjacent regions. To resolve this ambiguity, we perform a last step where we go through the different regions sequentially and assign boundary points to the current region unless previously assigned.

Figure 2 shows the obtained subdivision with $$K=18$$

, where it can be seen that the main regions of high curvature now belong to a unique subpart.

Fig. 2

The obtained Voronoi decomposition with 18 subparts

2.2 Conditional Model Parametrization

In the previous section we obtained $$K$$

subcomponents $\mathbf{x}_k$ , $k=1, \dots , K$ . The aim of this section is to describe the statistical modeling of the inter-part probability distributions, i.e. $P(\mathbf{x}_k|\mathbf{x}_l)$ , where $k, l=1, \dots , K$ and $k \ne l$ . More specifically, we would like to obtain new constraints for each part $\mathbf{x}_k$ based on its conditional relationship with $\mathbf{x}_l$ , that is, a new mean and covariance in the space of the shape parameters $\mathbf{b}_k$ . Let us denote ${\mu }_{kl}$ and ${\varSigma }_{kl}$ the values that form the new conditional constraints. In this paper, we choose to model $P(\mathbf{x}_k, \mathbf{x}_l)$ using a normal probability distribution. Thus, the mean and the covariance estimates are calculated as:

$\begin{aligned} {\mu }_{kl}={\varSigma }_{kl} {\varSigma }_{ll}^{-1}\mathbf {b}_l \end{aligned}$