Profiles: Simple, Fast, and Reliable Spine Localization in CT Scans

coordinate and pixel spacing, both given in millimeters. We denote the pixel positions in the patient coordinates as $\mathbf {p} = (p^x,p^y,p^z)$ . For sake of simplicity we also assume the feet-to-head, face-up (supine) orientation of the patient. The pixel gray values will be denoted as $g = g (\mathbf {p})$ and we assume them to be in Hounsfield units.

Our approach is based on bone distribution signatures within slices. For a fixed slice at location

we are therefore first interested in a rough segmentation $\fancyscript{B}_{z}$ of the bones. Given the Hounsfield intensities

, this can be achieved by an interval threshold, i.e., using two constants:

$\begin{aligned} \fancyscript{B}_{z} = \{\mathbf {p} = (p^x,p^y,p^z) \quad | \quad p^z = z \ \wedge \ g(\mathbf {p}) \in [400, 1050] \} \end{aligned}$

(1)

3.1 Centroids and Deviations

The simplest features are based on the centroid of segmentations $\fancyscript{B}_{z}$

$\begin{aligned} \mathbf {\mu _z} \ = \ \frac{1}{|\fancyscript{B}_{z}|} \sum _{\mathbf {p} \in \fancyscript{B}_{z}} \mathbf {p} \ = \ (\mu _z^x, \mu _z^y, z) \end{aligned}$

(2)

and on the length $\sigma _z$ of the associated standard deviation vector

$\begin{aligned} \sigma _z = \sqrt{ \frac{1}{|\fancyscript{B}_{z}|-1} \sum _{\mathbf {p} \in \fancyscript{B}_{z}} (p^x - \mu _z^x)^2 + (p^y - \mu _z^y)^2}. \end{aligned}$

(3)

Centroids $\mathbf {\mu _z}$ correlate with the spine reliably in the lumbar slices where pelvis, ribs, or head do not contribute to it.

The lumbar part can be characterized by deviation lengths $\sigma _z$ related to size of a vertebra seen in an axial slice [9]. Values of $\sigma _z$ larger than 40 mm indicate presence of non-vertebra bones.

Centroid refinement. While reliable in the lumbar area, the centroids $\mathbf {\mu _z}$ may drift remarkably from the spine if the pelvis or ribs contribute by their pixels (cf. Fig. 1).

Fig. 1

Examples showing centroids $\mathbf {\mu _z}$ (square), a circle of radius $\sigma _z$ , the $80 \times 140$ mm refinement window (rectangle), and refined center $\mathbf {\nu _z}$ (star)

To avoid this we refine the centroids $\mathbf {\mu _z}$ within rectangular $80 \times 140$ mm windows, asymmetrically spanned around them:

$\begin{aligned} \fancyscript{W}_{z} = \{\mathbf {p} \in \fancyscript{B}_{z} \,\, | \,\, -40 \le p^x - \mu _z^x \le 40 \,\,\, \wedge \,\, -40 \le p^y - \mu _z^y \le 100 \} \end{aligned}$

(4)

The size of the windows is set to be sufficiently big to accommodate any vertebra in an axial view [9] and to account for relative positions of the centroids and vertebrae in pelvis slices. The centroids $\mathbf {\mu _z}$ are refined to the center $\mathbf {\nu _z}$ of bone pixels in this window:

$\begin{aligned} \mathbf {\nu _z} = \frac{1}{|\fancyscript{W}_{z}|} \sum _{\mathbf {p} \in \fancyscript{W}_{z}} \mathbf {p} \end{aligned}$

(5)

3.2 Shape Histograms: AP Versus LR Distribution

To identify leg slices, we propose to discriminate slices with bone distributions dominant in the left-to-right direction and zero contributions in the anterior-posterior direction (see Fig. 2).

Fig. 2

Right, left, ante, and poste histogram bins centered at $\mathbf {\nu _z}$ overlaid over a negative of a CT slice at position $$z$$

showing legs and table

Fig. 3

An example CT scan in a frontal maximum intensity projection (MIP) (a), its histogram profile (b) and deviation profile (c)

We construct 4-bin histograms located in the refined centers $\mathbf {\nu _z}$ . Putting $\delta = \mathbf {p} - \mathbf {\nu _z}$ we define the following four quantities:

$\begin{aligned} h_z^A&= | \{ \mathbf {p} \in \fancyscript{B}_{z} \ | \ \delta ^y < - | \delta ^x | \le 0 \} | \end{aligned}$

(6)

$$$\begin{aligned} h_z^P&= | \{ \mathbf {p} \in \fancyscript{B}_{z} \ | \ \delta ^y > | \delta ^x | \ge 0 \} | \end{aligned}$$” src=”/wp-content/uploads/2016/10/A331518_1_En_15_Chapter_Equ7.gif”></DIV></DIV><br /> <DIV class=EquationNumber>(7)</DIV></DIV><br /> <DIV id=Equ8 class=Equation><br /> <DIV class=EquationContent><br /> <DIV class=MediaObject><IMG alt=$

(8)

$$$\begin{aligned} h_z^L&= | \{ \mathbf {p} \in \fancyscript{B}_{z} \ | \ \delta ^x > | \delta ^y | > 0 \} | \end{aligned}$$” src=”/wp-content/uploads/2016/10/A331518_1_En_15_Chapter_Equ9.gif”></DIV></DIV><br /> <DIV class=EquationNumber>(9)</DIV></DIV>With the AP/LR histograms we reformulate the leg detection as a search for slices, where ante-poste bone contributions vanish. We introduce scalars <SPAN id=IEq25 class=InlineEquation><IMG alt=$ and a threshold to yield this:

$\begin{aligned} 0 \le \lambda _z = \frac{h_z^A}{h_z^L + h_z^R} < 0.04 \end{aligned}$

(10)

Note that we exclude the posterior voxels ${h_z^P}$ from (10) in order to ignore eventual contribution of a CT table.

3.3 Bone Profiles

In the previous section we have introduced two bone distribution descriptors, i.e. scalars $\lambda _z$ and $\sigma _z$ for every slice $$z$$

. Next we aggregate them into two 1D arrays indexed by $$z$$

and refer to as the bone profiles. A symmetric plot of bone profiles along the z-axis is shown in Fig. 3.

4 Applications to Spinal Column Localization

In this section we show how the bone profiles and the refined centers $\mathbf {\nu _z}$ can be used to bound the spinal column and to identify a reliable initialization seed for subsequent computations.

4.1 Discarding the Slices up to the Ischium

When dealing with spines, leg slices should be taken out of consideration. We observed that the first occurrence of vanishing $\lambda _z$ in the top-to-bottom order may correspond either to the bottom of sacrum or the bottom of the pelvis—the ischium (cf. Fig. 3b).

In order to have a security margin between the spine and the slices to drop we suggest to identify ischium slices. We identify them by the first $$65$$

mm long segment of zeros in the histogram profile, i.e. a sequence longer than the average distance from ischium to bottom of sacrum.

4.2 Seeding a Spinal Canal Search

Algorithms using incremental/propagated search need to be initialized [5]. To obtain a reliable seed point near the spinal canal we consider the refined center $\mathbf {\nu _{z^\star }}$ in a slice with minimal deviation $\sigma _{z}$ (cf. Eq. 11). In this case no other bones except for vertebra contribute to the signatures and the point $\mathbf {\nu _{z^\star }}$ yields an estimate of the spinal canal. Such slices are predominantly found either in the lumbar area between pelvis and the first rib (cf. Fig. 3c) or in the neck area.

$\begin{aligned} \mathbf {\nu _{z^\star }} \quad | \quad z^\star = {\arg \!\min } \{ \sigma _z \} \end{aligned}$

(11)

4.3 Bounding the Spinal Column

Machine learning methods need to compute a vector of features at every voxel. Reducing the amount of voxels to be classified to a minimum can therefore significantly speed up such algorithms. After the leg slices have been discarded we wish to further prune the space by setting coronal and sagittal bounding planes.

For healthy spines the previously found seed $\mathbf {\nu _{z^\star }}$ could be reused to set up a bounding box of a predefined size. Such an approach would, however, fail for scolioses and other spine curvature related disorders.

To deal with such cases we derive coronal and sagittal planes from the bounding box of the $80 \times 140$ mm windows (cf. Eq. 4) spanned symmetrically around a subset of centers $\mathbf {\nu _z}$ . The refined centers $\mathbf {\nu _z}$ are first sorted by a drift reliability $\varDelta _{z}$ = $||\mathbf {\nu _z} - \mathbf {\mu _z}||$ from the original centroids $\mathbf {\mu _z}$ : the smaller the drift the more reliable the center. A fraction of sorted $\mathbf {\nu _z}$ involved in spanning bounding planes balances the tightness of bounding around spinal column and the data reduction. It is the last and the only free parameter in our method.