Determination of the Spine-Based Coordinate System for an Efficient Cross-Sectional Visualization of 3D Spine Images

of mutually orthogonal axes x, y and z, represented by unit vectors $${\hat{\mathbf{e}}}_{Ix} = [1,0,0]_{I}$$, $${\hat{\mathbf{e}}}_{Iy} = [0,1,0]_{I}$$ and $${\hat{\mathbf{e}}}_{Iz} = [0,0,1]_{I}$$, respectively. In the case of 3D spine images, the image-based coordinate system usually corresponds to the coordinate system of the imaging device (i.e. CT or MR scanner), in which the body is longitudinally aligned. The image-based coordinate system is, in general, aligned with the coordinate system of the body, and therefore its axes represent the following anatomical directions:




  • the sinistro–dexter axis x represents the direction $${\hat{\mathbf{e}}}_{Ix}$$ from the left to the right part of the body,


  • the ventro–dorsal axis y represents the direction $${\hat{\mathbf{e}}}_{Iy}$$ from the anterior to the posterior part of the body,


  • the cranio–caudal axis z represents the direction $${\hat{\mathbf{e}}}_{Iz}$$ from the superior to the inferior part of the body (i.e. the longitudinal axis of the body).


As a result, the following three imaging planes can be defined:



  • a sagittal or lateral plane $$(y,z)$$ is any plane that passes from the anterior to the posterior, and from the superior to the inferior part of the body, therefore dividing the body into left and right sections (the sagittal plane that passes through the midline of the body and, assuming bilateral symmetry, divides the body into left and right halves is the mid-sagittal plane),


  • a coronal or frontal plane $$(x,z)$$ is any plane that passes from the left to the right, and from the superior to the inferior part of the body, therefore dividing the body into anterior and posterior sections (the coronal plane that passes through the midline of the body and, assuming bilateral symmetry, divides the body into anterior and posterior halves is the mid-coronal plane),


  • an axial or transverse plane $$(x,y)$$ is any plane that passes from the left to the right, and from the anterior to the posterior part of the body, therefore dividing the body into superior and inferior sections.

The image-based coordinate system is therefore a right-handed Cartesian coordinate system with a standard orthonormal basis (Fig. 1).

A312884_1_En_8_Fig1_HTML.gif


Fig. 1
The image-based coordinate system $${\mathbb{R}}_{I}^{3} \to (x,y,z)$$ of a 3D image of a scoliotic spine, shown in a 3D view, b left sagittal view, c posterior coronal view and d superior axial view



2.1.2 Spine-Based Coordinate System


The spine-based coordinate system [84] is the coordinate system of the spine as the observed anatomical structure. In the case of 3D spine images, it is defined by a 3-tuple $$(u,v,w) \in {\mathbb{R}}_{S}^{3}$$ of mutually orthogonal axes $$u$$, $$v$$ and $$w$$, represented by unit vectors $${\hat{\mathbf{e}}}_{Su} = [1,0,0]_{S}$$, $${\hat{\mathbf{e}}}_{Sv} = [0,1,0]_{S}$$ and $${\hat{\mathbf{e}}}_{Sw} = [0,0,1]_{S}$$, respectively. The spine-based coordinate system is aligned with the spine, and its axes represent the following structural directions:



  • the anatomical axis $$u$$ represents the direction $${\hat{\mathbf{e}}}_{Su}$$ from the left to the right part of the spine,


  • the anatomical axis $$v$$ represents the direction $${\hat{\mathbf{e}}}_{Sv}$$ from the anterior to the posterior part of the spine,


  • the anatomical axis $$w$$ represents the direction $${\hat{\mathbf{e}}}_{Sw}$$ from the superior to the inferior part of the spine (i.e. the longitudinal axis of the spine).

As a result, the following three anatomical (structural) planes can be defined:



  • a sagittal or lateral plane $$(v,w)$$ is any plane that passes from the anterior to the posterior, and from the superior to the inferior part of the spine, therefore dividing the spine into left and right sections,


  • a coronal or frontal plane $$(u,w)$$ is any plane that passes from the left to the right, and from the superior to the inferior part of the spine, therefore dividing the spine into anterior and posterior sections,


  • an axial or transverse plane $$(u,v)$$ is any plane that passes from the left to the right, and from the anterior to the posterior part of the spine, therefore dividing the spine into superior and inferior sections.

The spine-based coordinate system is therefore a right-handed Cartesian coordinate system with a standard orthonormal basis (Fig. 2).

A312884_1_En_8_Fig2_HTML.gif


Fig. 2
The spine-based coordinate system $${\mathbb{R}}_{S}^{3} \to (u,v,w)$$ of a 3D image of a scoliotic spine, shown in a 3D view, b left sagittal view, c posterior coronal view and d superior axial view (Note The spine corresponds to Fig. 1)



2.2 Definition of the Spine-Based Coordinate System


To define the spine-based coordinate system, two geometrical properties of the spine have to be determined, namely the spine curve (Sect. 2.2.1) and axial vertebral rotation (Sect. 2.2.2), which serve to define the transformation from the image-based to the spine-based coordinate system (Sect. 2.2.3).


2.2.1 Spine Curve


The spine curve $$C$$ is the curve that follows the curvature of the spine along its entire longitudinal length. If $$i$$ is an independent parameter that denotes an arbitrary location on the spine, then $${\mathbf{c}}(i)$$ is the parametrization of the spine curve $$C$$:


$$C{:} \quad {\mathbf{c}}(i) = (c_{x} (i),c_{y} (i),c_{z} (i));\quad i \in [i_{\text{sp}} ,i_{\text{ep}} ],$$

(1)
where $$i = i_{\text{sp}}$$ and $$i = i_{\text{ep}}$$ represent the locations on the spine at its start and end point of observation, respectively, and $$c_{x} (i)$$, $$c_{y} (i)$$ and $$c_{z} (i)$$ represent the sagittal, coronal and axial coordinate, respectively, of the same anatomical reference point at any location $$i$$ on the spine in the image-based coordinate system. Although arbitrary anatomical reference points can be chosen (e.g. the centers of the spinal canal), the most established anatomical reference points are the centers of vertebral bodies. For $$K$$ observed consecutive vertebrae, let points $$\{ {\mathbf{v}}(k) = (v_{x} (k),v_{y} (k),v_{z} (k)); k = 1,2, \ldots ,K\}$$ represent the corresponding centers of vertebral bodies. The spine curve $${\mathbf{c}}(i)$$ can be then obtained by continuous interpolation of $${\mathbf{v}}(k)$$ between $${\mathbf{c}}(i_{\text{sp}} ) = {\mathbf{v}}(1)$$ and $${\mathbf{c}}(i_{\text{ep}} ) = {\mathbf{v}}(K)$$ (Fig. 3).

A312884_1_En_8_Fig3_HTML.gif


Fig. 3
The sagittal $$c_{x} (i)$$, coronal $$c_{y} (i)$$ and axial $$c_{z} (i)$$ component of the spine curve $${\mathbf{c}}(i)$$ against the independent parameter $$i \in [0,1]$$. Labels C7, T1, …, L3 indicate vertebral segments (Note The spine corresponds to Fig. 1)

For a differentiable curve $${\mathbf{c}}(i)$$, the geometric properties of the curve can be described in terms of differential geometry by the Frenet-Serret frame, which is defined as an orthonormal basis by the unit tangent, normal and binormal vectors to the curve. The unit tangent vector $${\hat{\mathbf{t}}}(i)$$ represents the direction of the curve that corresponds to increasing values of parameter $$i$$:


$${\hat{\mathbf{t}}}(i) = \frac{{{\mathbf{t}}(i)}}{{\left\| {{\mathbf{t}}(i)} \right\|}};\quad {\mathbf{t}}(i) = \frac{{{\text{d}}{\mathbf{c}}(i)}}{{{\text{d}}i}},$$

(2)
where $${\mathbf{t}}(i)$$ is the tangent vector to the curve $${\mathbf{c}}(i)$$, obtained as the first derivative of $${\mathbf{c}}(i)$$ with respect to $$i$$, and $$\left\| \cdot \right\|$$ denotes the vector norm. The unit normal vector $${\hat{\mathbf{n}}}(i)$$ represents the deviation of the curve from being a straight line:


$${\hat{\mathbf{n}}}(i) = \frac{{{\mathbf{n}}(i)}}{{\left\| {{\mathbf{n}}(i)} \right\|}};\quad {\mathbf{n}}(i) = \frac{{{\text{d}}{\hat{\mathbf{t}}}(i)}}{{{\text{d}}i}} = \frac{{{\text{d}}{\mathbf{c}}(i)}}{{{\text{d}}i}} \times \left( {\frac{{ {\text{d}^{2}}{\mathbf{c}}(i)}}{{{\text{d}}i^{2} }} \times \frac{{{\text{d}}{\mathbf{c}}(i)}}{{{\text{d}}i}}} \right),$$

(3)
where $${\mathbf{n}}(i)$$ is the normal vector to the curve $${\mathbf{c}}(i)$$, obtained as the first derivative of $${\hat{\mathbf{t}}}(i)$$ with respect to $$i$$, and $$\times$$ denotes the cross vector product. To satisfy the orthonormality of the basis, the unit binormal vector $${\hat{\mathbf{b}}}(i)$$ is orthogonal to both the unit tangent vector and the unit normal vector:


$${\hat{\mathbf{b}}}(i) = {\hat{\mathbf{t}}}(i) \times {\hat{\mathbf{n}}}(i) = \frac{{{\mathbf{b}}(i)}}{{\left\| {{\mathbf{b}}(i)} \right\|}};\quad {\mathbf{b}}(i) = \frac{{{\text{d}}{\mathbf{c}}(i)}}{{{\text{d}}i}} \times \frac{{{\text{d}}^{2} {\mathbf{c}}(i)}}{{{\text{d}}i^{2} }},$$

(4)
where $${\mathbf{b}}(i)$$ is the binormal vector to the curve $${\mathbf{c}}(i)$$, obtained as the cross vector product of the first and the second derivative of the curve $${\mathbf{c}}(i)$$.

The unit tangent vector $${\hat{\mathbf{t}}}(i)$$ and the unit normal vector $${\hat{\mathbf{n}}}(i)$$ at location $$i$$ on curve $${\mathbf{c}}(i)$$ define the osculating plane at that location. The deviation of the curve from being a straight line relative to the osculating plane is measured by the geometrical curvature $${\kappa} (i)$$ (Fig. 4):

A312884_1_En_8_Fig4_HTML.gif


Fig. 4
The geometrical curvature $${\kappa}(i)$$ of the spine curve $${\mathbf{c}}(i)$$ against the independent parameter $$i \in [0,1]$$. Labels C7, T1, …, L3 indicate vertebral segments (Note The spine corresponds to Fig. 1)



$${\kappa} (i) = \frac{1}{{r_{\kappa} (i)}} = \frac{{\left\| {\frac{{{\text{d}}{\mathbf{c}}(i)}}{{{\text{d}}i}} \times \frac{{{\text{d}}^{2} {\mathbf{c}}(i)}}{{{\text{d}}i^{2} }}} \right\|}}{{\left\| {\frac{{{\text{d}}{\mathbf{c}}(i)}}{{{\text{d}}i}}} \right\|^{3} }},$$

(5)
where $$r_{\kappa } (i)$$ is the radius of curvature that represents the radius of the osculating circle in the osculating plane. On the other hand, the deviation of the curve from being a plane curve, represented by the rotation of the unit binormal vector $${\hat{\mathbf{b}}}(i)$$ about the unit tangent vector $${\hat{\mathbf{t}}}(i)$$, is measured by the geometrical torsion $$\tau (i)$$ (Fig. 5):

A312884_1_En_8_Fig5_HTML.gif


Fig. 5
The geometrical torsion $$\tau (i)$$ of the spine curve $${\mathbf{c}}(i)$$ against the independent parameter $$i \in [0,1]$$. Labels C7, T1, …, L3 indicate vertebral segments (Note The spine corresponds to Fig. 1)



$$\tau (i) = \frac{1}{{r_{\sigma } (i)}} = \frac{{\left( {\frac{{{\text{d}}{\mathbf{c}}(i)}}{{{\text{d}}i}} \times \frac{{{\text{d}}^{2} {\mathbf{c}}(i)}}{{{\text{d}}i^{2} }}} \right) \cdot \frac{{{\text{d}}^{3} {\mathbf{c}}(i)}}{{{\text{d}}i^{3} }}}}{{\left\| {\frac{{{\text{d}}{\mathbf{c}}(i)}}{{{\text{d}}i}} \times \frac{{{\text{d}}^{2} {\mathbf{c}}(i)}}{{{\text{d}}i^{2} }}} \right\|^{2} }},$$

(6)
where $$\cdot$$ denotes the dot vector product, and $$r_{\sigma } (i)$$ is the radius of torsion. By using the geometrical curvature and torsion, the resulting Frenet-Serret frame can be written in matrix form as:


$$\frac{\text{d}}{{{\text{d}}i}}\left[ {\begin{array}{*{20}c} {{\hat{\mathbf{t}}}(i)} \\ {{\hat{\mathbf{n}}}(i)} \\ {{\hat{\mathbf{b}}}(i)} \\ \end{array} } \right] = \left\| {\frac{{{\text{d}}{\mathbf{c}}(i)}}{{{\text{d}}i}}} \right\|\left[ {\begin{array}{*{20}c} 0 & {\kappa (i)} & 0 \\ { - \kappa (i) } & 0 & { \tau (i) } \\ 0 & { - \tau (i) } & 0 \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {{\hat{\mathbf{t}}}(i)} \\ {{\hat{\mathbf{n}}}(i)} \\ {{\hat{\mathbf{b}}}(i)} \\ \end{array} } \right].$$

(7)

The form of Eqs. 2–7 corresponds to a regular parametrization of the curve by location i on the curve (Eq. 1). In the case the curve is reparameterized by its arc length s, the natural parametrization $${\mathbf{c}}(s)$$ of $${\mathbf{c}}(i)$$ is yielded:


$${\mathbf{c}}(s) = {\mathbf{c}}(i(s));\quad i(s) = s^{ - 1} (i);\quad s(i) = \int\limits_{{i_{\text{sp}} }}^{i} \left\| {\frac{{{\text{d}}{\mathbf{c}}(\lambda )}}{{{\text{d}}\lambda }}} \right\|{\text{d}}\lambda .$$

(8)

Considering the natural parametrization of the curve, the unit tangent vector $${\hat{\mathbf{t}}}(s)$$, unit normal vector $${\hat{\mathbf{n}}}(s)$$ and unit binormal vector $${\hat{\mathbf{b}}}(s)$$ are computed as:


$${\hat{\mathbf{t}}}(s) = \frac{{{\text{d}}{\mathbf{c}}(s)}}{{{\text{d}}s}}; \quad {\hat{\mathbf{n}}}(s) = \frac{{\frac{{{\text{d}\hat{\mathbf{t}}}(s)}}{{{\text{d}}s}}}}{{\left\| {\frac{{{\text{d}}{\hat{\mathbf{t}}}(s)}}{{{\text{d}}s}}} \right\|}}; \quad {\hat{\mathbf{b}}}(s) = {\hat{\mathbf{t}}}(s) \times {\hat{\mathbf{n}}}(s),$$

(9)
the corresponding geometrical curvature $${\kappa} (s)$$ and torsion $$\tau (s)$$ are computed as:


$$\kappa (s) = \left\| {\frac{{{\text{d}}{\hat{\mathbf{t}}}(s)}}{{{\text{d}}s}}} \right\| = \left\| {\frac{{{\text{d}}^{2} {\mathbf{c}}(s)}}{{{\text{d}}s^{2} }}} \right\|; \quad \tau (s) = - {\hat{\mathbf{n}}{(\mathbf s)}} \cdot \frac{{{\text{d}}{\hat{\mathbf{b}}}(s)}}{{{\text{d}}s}},$$

(10)
and the Frenet-Serret frame in the matrix form is:


$$\frac{\text{d}}{\text{d}{s}}\left[ {\begin{array}{*{20}c} {{\hat{\mathbf{t}}}(s)} \\ {{\hat{\mathbf{n}}}(s)} \\ {{\hat{\mathbf{b}}}(s)} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} 0 & {\kappa (s)} & 0 \\ { - \kappa (s) } & 0 & { \tau (s) } \\ 0 & { - \tau (s) } & 0 \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {{\hat{\mathbf{t}}}(s)} \\ {{\hat{\mathbf{n}}}(s)} \\ {{\hat{\mathbf{b}}}(s)} \\ \end{array} } \right].$$

(11)

However, the natural parametrization is in the case of spine curves rare, as it is often based on the axial coordinate $$z$$ at the start and end point of the spine curve (i.e. $$i_{\text{sp}} = z_{1}$$ and $$i_{\text{ep}} = z_{2}$$, where $${\mathbf{p}}_{1} = (x_{1} ,y_{1} ,z_{1} )$$ is the start point and $${\mathbf{p}}_{2} = (x_{2} ,y_{2} ,z_{2} )$$ is the end point), or on pre-defined constant values (e.g. $$i_{\text{sp}} = 0$$ and $$i_{\text{ep}} = 1$$).


2.2.2 Axial Vertebral Rotation


Axial vertebral rotation is the rotation $${\varPhi }$$ of vertebrae around their longitudinal axes when projected onto the transverse plane of the corresponding coordinate system. If $$i$$ is an independent parameter that denotes an arbitrary location on the spine, then $$\varphi (i)$$ is the parametrization of the axial vertebral rotation $${\varPhi }$$:


$${\varPhi } {:} \quad \varphi (i);\quad i \in [i_{\text{sp}} ,i_{\text{ep}} ],$$

(12)
where $$i = i_{\text{sp}}$$ and $$i = i_{\text{ep}}$$ represent the locations on the spine at its start and end point of observation, respectively. For $$K$$ observed consecutive vertebrae, let angles $$\{ \delta (k); k = 1,2, \ldots ,K\}$$ represent the corresponding axial vertebral rotation angles. The axial vertebral rotation $$\varphi (i)$$ can be then obtained by continuous interpolation of $$\delta (k)$$ between $$\varphi (i_{\text{sp}} ) = \delta (1)$$ and $$\varphi (i_{\text{ep}} ) = \delta (K)$$ (Fig. 6).

A312884_1_En_8_Fig6_HTML.gif


Fig. 6
The axial vertebral rotation $$\varphi (i)$$ against the independent parameter $$i \in [0,1]$$. Labels C7, T1, …, L3 indicate vertebral segments (Note The spine corresponds to Fig. 1)

Axial vertebral rotation cannot be uniformly measured, as the vertebral anatomy observed in the transverse plane is not completely symmetrical due to normal developmental as well as pathological conditions affecting the spine. As a result, several methods were developed to determine axial vertebral rotation from 3D images of the spine [88] that measured the angle between the reference sagittal plane and a line connecting specific anatomical reference points in the transverse plane. For example, axial vertebral rotation was measured as the angle between the reference sagittal plane and the line connecting the posterior junction of the two laminae of the vertebral arch with the center of vertebral body [1], as the angle between the reference sagittal plane and the line bisecting the angle between the two lines connecting the junction of each lamina and the pedicle with the posterior junction of the two laminae [26], as the angle between the reference sagittal plane and the line connecting the tip of the spinous process with the center of vertebral body [39], or as the angle between the reference sagittal plane and the line connecting the most posterior points of the two pedicles [20].

Let points $${\mathbf{r}}(i) = (r_{x} (i),r_{y} (i),r_{z} (i))$$ represent locations of the selected anatomical reference points in the corresponding transverse planes of measurement, and let $${\mathbf{c}}(i) = (c_{x} (i),c_{y} (i),c_{z} (i))$$ represent the spine curve that passes through the centers of vertebral bodies. If the transverse planes of measurement are image-based, i.e. orthogonal to axis $$z$$ of the image-based coordinate system, then the axial vertebral rotation $$\varphi (i) = \varphi_{z} (i)$$ can be determined by considering $$r_{z} (i) = c_{z} (i)$$ as:


$$\varphi_{z} (i) = \arctan \left( {\frac{{r_{x} (i) - c_{x} (i)}}{{r_{y} (i) - c_{y} (i)}}} \right),$$

(13)
which for every $$i$$ corresponds to the angle between the line connecting the anatomical reference point with the center of the vertebral body, and the line representing the reference sagittal plane (i.e. the line in the direction of axis $$y$$ of the image-based coordinate system). However, because vertebrae can be sagittally or coronally inclined against axis $$z$$, the centers of vertebral bodies and the corresponding reference anatomical points may not represent corresponding anatomical locations along the longitudinal vertebral axes. On the other hand, if the transverse planes of measurement are spine-based, i.e. orthogonal to axis $$w$$ of the spine-based coordinate system and therefore orthogonal to $${\mathbf{c}}(i)$$, then the axial vertebral rotation $$\varphi (i) = \varphi_{w} (i)$$ is measured at corresponding anatomical locations and can be determined as:


$$\varphi_{w} (i) = \arccos \left( {\frac{{({\mathbf{r}}(i) - {\mathbf{c}}(i)) \cdot {\tilde{\mathbf{e}}}_{Iy} (i)}}{{\left\| {{\mathbf{r}}(i) - {\mathbf{c}}(i)} \right\|}}} \right),\quad {\tilde{\mathbf{e}}}_{Iy} (i) = {\hat{\mathbf{t}}}(i) \times \left( {{\hat{\mathbf{e}}}_{Iy} \times {\hat{\mathbf{t}}}(i)} \right),$$

(14)
where $${\tilde{\mathbf{e}}}_{Iy} (i)$$ is the unit vector in the direction of the projection of $${\hat{\mathbf{e}}}_{Iy} = [0,1,0]_{I}$$ to the plane orthogonal to the spine curve, defined by the unit tangent vector $${\hat{\mathbf{t}}}(i)$$ as the normal of that plane.


2.2.3 Transformation from Image-Based to Spine-Based Coordinate System


The continuous transformation from the image-based coordinate system (Fig. 7) to the spine-based coordinate system (Fig. 8) is possible by having a continuous description of the spine curve $$C$$ (Eq. 1) and axial vertebral rotation $${\varPhi }$$ (Eq. 12) that are parameterized by the same variable $$i$$ representing the location on the spine:

A312884_1_En_8_Fig7_HTML.gif


Fig. 7
The spine curve (red line) and axial vertebral rotation (blue directions) in the image-based coordinate system $${\mathbb{R}}_{I}^{3}$$ of a 3D image of a scoliotic spine, shown in a 3D view, b left sagittal view, c posterior coronal view and d superior axial view. The spine curve is defined between the start point $${\mathbf{p}}_{1}$$ and end point $${\mathbf{p}}_{2}$$, while the coordinates are shown also for a selected point $${\mathbf{p}}_{c} = (x_{c} ,y_{c} ,z_{c} )$$ on the spine curve


A312884_1_En_8_Fig8_HTML.gif


Fig. 8
The spine curve (red line) and axial vertebral rotation (blue directions) in the spine-based coordinate system $${\mathbb{R}}_{S}^{3}$$ of a 3D image of a scoliotic spine, shown in a 3D view, b left sagittal view, c posterior coronal view and d superior axial view. The spine curve is defined between the start point $${\mathbf{p}}_{1}$$ and end point $${\mathbf{p}}_{2}$$, while the coordinates are shown also for a selected point $${\mathbf{p}}_{c} = (u_{c} ,v_{c} ,w_{c} )$$ on the spine curve



$$\{ C, {\varPhi } \} {:} \quad \zeta (i) = ({\mathbf{c}}(i),\varphi (i)) = (c_{x} (i),c_{y} (i),c_{z} (i),\varphi (i));\quad i = [i_{\text{sp}} ,i_{\text{ep}} ],$$

(15)
where $$i = i_{\text{sp}}$$ and $$i = i_{\text{ep}}$$ represent the locations of the spine curve and axial vertebral rotation at the start and end point of observation on the spine, respectively.

The Frenet-Serret frame (Eq. 7) describes the geometrical properties of the curve and can be therefore applied to the spine curve $${\mathbf{c}}(i)$$. However, the spine-based coordinate system $$(u,v,w) \in {\mathbb{R}}_{S}^{3}$$ has to represent also the course of the axial vertebral rotation $$\varphi (i)$$. The unit tangent vector $${\hat{\mathbf{t}}}(i)$$ defines the unit vector $${\hat{\mathbf{e}}}_{Sw} = [0,0,1]_{S}$$ representing axis $$w$$ at any point $$i$$ on the spine curve. At the same time it defines, as the normal, the plane orthogonal to the spine curve at any point $$i$$, which is the plane that contains unit vectors $${\hat{\mathbf{e}}}_{Su} = [1,0,0]_{S}$$ and $${\hat{\mathbf{e}}}_{Sv} = [0,1,0]_{S}$$ representing axes $$u$$ and $$v$$, respectively. However, $${\hat{\mathbf{e}}}_{Su}$$ and $${\hat{\mathbf{e}}}_{Sv}$$ do not correspond to the unit normal vector $${\hat{\mathbf{n}}}(i)$$ or unit binormal vector $${\hat{\mathbf{b}}}(i)$$, although they lie in the same plane, i.e. the plane orthogonal to the spine curve, because $${\hat{\mathbf{n}}}(i)$$ and $${\hat{\mathbf{b}}}(i)$$ rotate about $${\hat{\mathbf{t}}}(i)$$ as described by the geometrical torsion $$\tau (i)$$. The directions of $${\hat{\mathbf{e}}}_{Su}$$ and $${\hat{\mathbf{e}}}_{Sv}$$ are namely defined by the axial vertebral rotation $$\varphi (i)$$. If $$\varphi (i) = \varphi_{z} (i)$$ (Eq. 13), meaning that the axial vertebral rotation is defined in transverse planes of measurement that are orthogonal to axis $$z$$ of the image-based coordinate system, then the modified unit normal vector $${\hat{\mathbf{n}}}_{\varphi } (i)$$ equals:


$${\hat{\mathbf{n}}}_{\varphi } (i) = {\hat{\mathbf{t}}}(i) \times (R_{z} (\varphi_{z} (i)) {\hat{\mathbf{e}}}_{Iy} \times {\hat{\mathbf{t}}}(i)),$$

(16)
where matrix $$R_{z} (\varphi_{z} (i))$$ represents the rotation for angle $$\varphi_{z} (i)$$ about axis $$z$$ of the image-based coordinate system:


$$R_{z} (\varphi_{z} (i)) = \left[ {\begin{array}{*{20}c} {\cos \varphi_{z} (i)} & { - \sin \varphi_{z} (i)} & 0 \\ {\sin \varphi_{z} (i)} & ~~~~{\cos \varphi_{z} (i)} & 0 \\ {\;\;\;\,0} & {\;\;\,0} & 1 \\ \end{array} } \right],$$

(17)
with the center of rotation at point $${\mathbf{c}}(i) = (c_{x} (i),c_{y} (i),c_{z} (i))$$. On the other hand, if $$\varphi (i) = \varphi_{w} (i)$$ (Eq. 14), meaning that the axial vertebral rotation is defined in transverse planes of measurement that are orthogonal to the spine curve $${\mathbf{c}}(i)$$, then the modified unit normal vector $${\hat{\mathbf{n}}}_{\varphi } (i)$$ equals:


$${\hat{\mathbf{n}}}_{\varphi } (i) = R_{{{\hat{\mathbf{t}}}(i)}} (\varphi_{w} (i)) {\tilde{\mathbf{e}}}_{Iy} (i),$$

(18)
where $${\tilde{\mathbf{e}}}_{Iy} (i)$$ is the unit vector in the direction of the projection of $${\hat{\mathbf{e}}}_{Iy} = [0,1,0]_{I}$$ to the plane orthogonal to the spine curve (Eq. 14), and matrix $$R_{{{\hat{\mathbf{t}}}(i)}} (\varphi_{w} (i))$$ represents the rotation for angle $$\varphi_{w} (i)$$ about the axis defined by the unit tangent vector $${\hat{\mathbf{t}}}(i)$$:


$$\begin{aligned} R_{{{\hat{\mathbf{t}}}(i)}} (\varphi_{w} (i)) & = \cos (\varphi_{w} (i)) I_{3} \\ & \quad + \sin (\varphi_{w} (i)) [{\hat{\mathbf{t}}}(i)]_{ \times } \\ & \quad + (1 - \cos (\varphi_{w} (i))) [{\hat{\mathbf{t}}}(i) \otimes {\hat{\mathbf{t}}}(i)], \\ \end{aligned}$$

(19)
with the center of rotation at point $${\mathbf{c}}(i) = (c_{x} (i),c_{y} (i),c_{z} (i))$$. In Eq. 19, $$I_{3}$$ denotes the identity matrix of size $$3{\kern 1pt} \times {\kern 1pt} 3$$, and $$[{\hat{\mathbf{t}}}(i)]_{ \times }$$ and $$[{\hat{\mathbf{t}}}(i){\kern 1pt} \otimes {\kern 1pt} {\hat{\mathbf{t}}}(i)]$$ are, respectively, the cross and tensor product matrix of $${\hat{\mathbf{t}}}(i)$$:


$$[{\hat{\mathbf{t}}}(i)]_{ \times } = \left[ {\begin{array}{*{20}c} ~~~0 & { - \hat{t}_{z} (i) } & ~~~{ \hat{t}_{y} (i) } \\ ~~~{ \hat{t}_{z} (i) } & 0 & { - \hat{t}_{x} (i) } \\ { - \hat{t}_{y} (i) } & ~~~{ \hat{t}_{x} (i) } & ~~~0 \\ \end{array} } \right],$$

(20)



$$[{\hat{\mathbf{t}}}(i) \otimes {\hat{\mathbf{t}}}(i)] = \left[ {\begin{array}{*{20}c} { (\hat{t}_{x} (i))^{2} } & { \hat{t}_{x} (i){\kern 1pt} \hat{t}_{y} (i) } & { \hat{t}_{x} (i){\kern 1pt} \hat{t}_{z} (i) } \\ { \hat{t}_{x} (i){\kern 1pt} \hat{t}_{y} (i) } & { (\hat{t}_{y} (i))^{2} } & { \hat{t}_{y} (i){\kern 1pt} \hat{t}_{z} (i) } \\ { \hat{t}_{x} (i){\kern 1pt} \hat{t}_{z} (i) } & { \hat{t}_{y} (i){\kern 1pt} \hat{t}_{z} (i) } & { (\hat{t}_{z} (i))^{2} } \\ \end{array} } \right].$$

(21)

In both cases, the unit binormal vector also changes its direction to fit the orthonormal basis, and is therefore equal to $${\hat{\mathbf{b}}}_{\varphi } (i) = {\hat{\mathbf{t}}}(i){\kern 1pt} \times {\kern 1pt} {\hat{\mathbf{n}}}_{\varphi } (i)$$. The resulting axes $$u$$, $$v$$ and $$w$$ of the spine-based coordinate system are therefore represented by:


$$u {:} \quad ~~~{\hat{\mathbf{e}}}_{Su} = [1,0,0]_{S} \quad \leftrightarrow \quad {\hat{\mathbf{b}}}_{\varphi } (i),$$

(22)



$$v {:} \quad ~~~{\hat{\mathbf{e}}}_{Sv} = [0,1,0]_{S} \quad \leftrightarrow \quad {\hat{\mathbf{n}}}_{\varphi } (i),$$

(23)



$$w {:} \quad {\hat{\mathbf{e}}}_{Sw} = [0,0,1]_{S} \quad \leftrightarrow \quad {\hat{\mathbf{t}}}(i).$$

(24)

In contrast to the image-based coordinate system, which is defined in the Euclidean space, the spine-based coordinate system is defined in a non-Euclidean space, and therefore morphometric measurements based on Euclidean metrics can not be obtained directly from the spine-based coordinate system.


2.3 Automated Determination of the Spine-Based Coordinate System


The spine-based coordinate system can be manually determined by identifying distinctive anatomical points on each vertebra (e.g. the centers of vertebral bodies) and the corresponding rotation of vertebrae, and then interpolating through these points to obtain a continuous description of both the spine curve and axial vertebral rotation along the whole length of the spine. However, navigation through 3D spine images is time-consuming and subjective, moreover it is practically impossible to manually define the plane orthogonal to the spine curve, in which axial vertebral rotation is defined, basing only on the identified anatomical points on each vertebra. As a result, several automated and semi-automated methods based on image processing and analysis techniques were proposed to determine the spine curve and/or axial vertebral rotation in 3D spine images.


2.3.1 Automated Determination of the Spine Curve


In the past, the only possibility for measuring the geometrical properties of the spine curve was based on examining the antero-posterior and/or lateral radiographs. As a result, the spine curve in 3D was observed as its projection in 2D in the form of coronal and sagittal spinal curvatures. Moreover, these curvatures were usually evaluated by one-dimensional measures including angles of curvature (e.g. the Ferguson angle, the Cobb angle, the tangent line angle) and indices of curvature (e.g. the Greenspan index, the Ishihara index). A detailed review of methods for the determination of spinal curvature was performed by Vrtovec et al. [89].

With the development of 3D imaging techniques, methods that captured the 3D nature of the spine started to emerge, followed by application of computerized techniques that automatically or semi-automatically (i.e. with minimal observer interaction) determined the spine curve in 3D images. Due to the continuous course of the spinal curvature, a number of studies attempted to model the spine curve with a mathematical curve in stereoradiographic (i.e. in multiple radiographs acquired at different angles), CT or MR spine images. Different functions were used for modeling, such as harmonic functions (i.e. sines, cosines or Fourier series) [13, 15, 27, 58, 74], spline functions [6, 33, 55, 82] and standard polynomial functions [56, 83, 85, 86], as well as statistical interpolation techniques, such as kriging [59].

By computerized least-squares aligning of a parametric sine function to the stereographically reconstructed landmarks, Stokes et al. [74] measured the Cobb angle between the normals to the obtained curve at inflection points in the coronal and sagittal plane, and in the plane of maximal curvature. Drerup and Hierholzer [13–15] also considered the sine function appropriate, as it most resembled the appearance of curves in idiopathic scoliosis. On the other hand, Patwardhan et al. [55] justified the use of spline functions by stating that splines are used to describe geometries with continuously changing curvature, such as scoliotic spines. In their framework for spine segmentation from CT images, Kaminsky et al. [33] used spline functions because they proved appropriate to describe both the anatomical shape and scoliotic deformations of the spine. Berthonnaud and Dimnet [6] constructed the spine curve separately in coronal and sagittal projections by computing the average of two spline functions that connected the anatomical landmarks on vertebral body walls. On the other hand, Peng et al. [56] used polynomial functions to detect and segment vertebrae from MR images using vertebral disc templates. Polynomial functions were also used to model both normal and pathological spine curves in CT images by Vrtovec et al. [83]. The spine curve was automatically determined by aligning the polynomial function with the centers of vertebral bodies in 3D, obtained by maximizing the distance from the edges of vertebral bodies. The same authors also developed a method for MR images [86], where the center of vertebral body was first automatically detected in each axial cross-section by maximizing the entropy of image intensities inside a circular region, and the detected centers of vertebral bodies in 3D were then joined by a polynomial function using the robust least-trimmed-squares regression. The method was also used by Neubert et al. [49] for extracting the spine curve from MR images of high resolution, obtained by applying the sequence named sampling perfection with application optimized contrasts using different flip angle evolution (SPACE). The work was continued by Štern et al. [78], who proposed a modality-independent method for the determination of the spine curve that was extracted from locations where lines connecting opposite edge points on vertebral body walls in the direction of corresponding image intensity gradients most often intersected. Kadoury et al. [30, 31] determined the spine curve of a scoliotic spine in biplanar radiographs by first embedding cubic B-spline functions onto a non-linear manifold to predict an initial curve according to a given database of scoliotic curves, and then performing analytical regression to obtain a statistical model of the final curve.

To extract the spinal canal centerline from CT images, Yao et al. [96] applied a watershed algorithm followed by a graph search, while Hay et al. [24] applied the fast marching minimal path technique that was based on the distance transform of the spinal canal segmentation, obtained by morphological region growing. On the other hand, Klinder et al. [36] segmented the spinal canal by a progressive adaptation of small tubular segments, represented as triangulated surface meshes, and then determined the spinal canal centerline by calculating the centers of mass for all contours of the obtained tubular mesh. A similar approach was proposed by Forsberg et al. [17, 18], who extracted the spinal canal centerline in CT images by first modeling the vertebral foramen in 2D axial cross-sections with circles, and then fitting a cubic B-spline function to the centers of the obtained circles.

Besides modeling the spine curve in 3D, different geometrical descriptors of spinal curvature were derived from mathematical functions. Poncet et al. [58] proposed the geometrical torsion as a measure for classifying spinal deformities, and Kadoury et al. [29, 32] further showed that it can be potentially used to discriminate among different types of thoracolumbar deformations of the spine. Vrtovec et al. [85] showed that clinically relevant features of the spine can be identified in 3D by observing the geometrical curvature as well as the curvature angle, which was defined as the angular magnitude of the geometrical curvature on an arbitrary spine section, and was as such independent of the size of the spine. Hay et al. [24] observed both the geometrical curvature and geometrical torsion that were scaled to a subject-independent coordinate system, and showed that they can be used to detect and quantify pathological spinal curvatures.


2.3.2 Automated Determination of the Axial Vertebral Rotation


Similarly as for the spinal curvature, measurements of axial vertebral rotation were in the past possible only by examining the location of pedicles and spinous processes in relation to corresponding vertebral bodies in antero-posterior radiographs. As a result, the axial vertebral rotation in 3D was observed as its projection in 2D, and several methods based on indices (e.g. the Cobb method, the Nash-Moe method, the Fait-Janovec method) or actual angles (e.g. the Bunnell method, the Drerup method, the Stokes et al. method) were proposed. With the introduction of 3D imaging techniques, cross-sectional imaging in the axial plane became possible and stimulated the development of methods that were based on manual identification of distinctive anatomical reference points (e.g. the tip of spinous process, the center of the vertebral body, etc.). A detailed review of methods for the determination of axial vertebral rotation was performed by Lam et al. [40] and Vrtovec et al. [88].

The measurement of axial vertebral rotation was also approached by computerized techniques based on image processing and analysis, although manual initialization was still required. Haughton et al. [23] and Rogers et al. [66, 67] proposed a method that required manual determination of the axial CT [66] or MR [23, 67] cross-section, the center of rotation and the circular area that encompassed the observed lumbar vertebra. After initialization, the method automatically measured the axial vertebral rotation relative to the cross-section of a different vertebra by searching for the maximal correlation of image intensities between the circular areas determined in both cross-sections. Oblique CT cross-sections were used by Adam and Askin [2], who determined the axial vertebral rotation from the line that bisected the thresholded image of the vertebral body according to the symmetry ratio, defined by the maximal correlation of image intensities in the bisected regions. Kouwenhoven et al. [37, 38] manually selected axial cross-sections through the centers of vertebral bodies in CT [38] and MR [37] images of normal spines, and applied automated region growing segmentation to obtain reference points, such as the center of the vertebral canal, the center of the sternum at the T5 vertebra and the center of the anterior half of the vertebral body, which were used to define the axial vertebral rotation. Axial vertebral rotation was studied in both CT and MR images of whole spines also by Vrtovec et al. [83, 86]. In CT images [83], circular cross-sections that were orthogonal to the spine curve were first automatically extracted, and axial vertebral rotation was then defined from the line that bisected the cross-section and resulted in the maximal correlation of image intensities in the bisected regions. For MR images [86], the rotation was defined in an optimization procedure that searched for the orientation angle of the line of symmetry in each axial-cross section, and then smoothed with a polynomial function along the whole spine using the least-trimmed-squares regression technique. The same authors also combined both approaches into a method that was modality-independent, i.e. applicable to both CT and MR images [87]. Basing on the pre-defined location of the vertebral body center in 3D, they obtained the relation between the image-based and vertebra-based coordinate systems by matching image intensity gradients that defined the best available symmetry of vertebral anatomical structures. The method was thoroughly evaluated and compared to established manual methods when applied to CT [91] and MR [90] images of normal and scoliotic spines. To segment vertebral bodies in both CT and MR images, Štern et al. [79] proposed to use a parametric model based on superquadrics that, among several shape parameters, included also the axial rotation of the vertebral body, and which was later used to perform quantitative vertebral morphometry in CT images of normal and fractured vertebrae [80]. Axial vertebral rotation was determined from the symmetry of vertebral anatomical structures also in the study of Forsberg et al. [18], who for each vertebra in CT images extracted a cross-section that was orthogonal to the spine curve and passed through the center of the vertebral body, and then minimized the sum of absolute differences in image intensities over the line that bisected the cross-section at the evaluated rotation angle. The same group of authors also developed a method for segmentation of vertebrae by registering a spine model to CT spine images, and then measured axial vertebral rotation from landmarks that were placed at distinctive anatomical locations in the spine model and mapped to each CT image by using the obtained registration transformation fields [17].


2.3.3 Examples of Automated Determination of the Spine Curve and Axial Vertebral Rotation


Among automated methods for the determination of the spine curve $${\mathbf{c}}(i)$$ and/or axial vertebral rotation $$\varphi (i)$$, the following approaches are presented in detail:


In all of the presented approaches [78, 83, 86], the spine-based coordinate system is determined from 3D spine images of normal and scoliotic subjects by parameterizing the spine curve $${\mathbf{c}}(i) = (c_{x} (i),c_{y} (i),c_{z} (i))$$ and axial vertebral rotation $$\varphi (i)$$ with polynomial functions:


$${\mathbf{c}}(i) = \left( {\sum\limits_{k = 0}^{{K_{x} }} \frac{{b_{x,k} }}{{\hat{b}_{k} }}{\kern 1pt} i^{{{\kern 1pt} k}} ,\sum\limits_{k = 0}^{{K_{y} }} \frac{{b_{y,k} }}{{\hat{b}_{k} }}{\kern 1pt} i^{{{\kern 1pt} k}} ,\sum\limits_{k = 0}^{{K_{z} }} \frac{{b_{z,k} }}{{\hat{b}_{k} }}{\kern 1pt} i^{{{\kern 1pt} k}} } \right),$$

(25)



$$\varphi (i) = \sum\limits_{k = 0}^{{K_{\varphi } }} \frac{{b_{{_{\varphi } ,k}} }}{{\hat{b}_{k} }}{\kern 1pt} i^{{{\kern 1pt} k}} ,$$

(26)
where $$b_{x} = \{ b_{x,k} ;k = 0,1, \ldots ,K_{x} \}$$, $$b_{y} = \{ b_{y,k} ;k = 0,1, \ldots ,K_{y} \}$$ and $$b_{z} = \{ b_{z,k} ;k = 0,1, \ldots ,K_{z} \}$$ are the parameters of polynomial functions $$c_{x} (i)$$, $$c_{y} (i)$$ and $$c_{z} (i)$$ of degrees $$K_{x}$$, $$K_{y}$$ and $$K_{z}$$, respectively, corresponding to the spine curve $${\mathbf{c}}(i)$$, and $$b_{\varphi } = \left\{ {b_{\varphi ,k} ;k = 0,1, \ldots ,K_{\varphi } } \right\}$$ are the parameters of the polynomial function $$\varphi (i)$$ of the degree $$K_{\varphi }$$ corresponding to the axial vertebral rotation $$\varphi (i)$$. The normalization coefficients $$\hat{b}_{k}$$:


$$\hat{b}_{k} = \int\limits_{{i_{\text{sp}} }}^{{i_{\text{ep}} }} |i^{{{\kern 1pt} k}} |{\text{d}}{i},$$

(27)
where $$i = i_{\text{sp}}$$ and $$i = i_{\text{ep}}$$ represent the locations on the spine at its start and end point of observation, respectively, regularize the impact of each polynomial parameter to the absolute variation of the corresponding term. With such parametrization, the goal is to automatically determine polynomial parameters $$b_{{\mathbf{c}}} = b_{x} \cup b_{y} \cup b_{z}$$ that describe the spine curve $${\mathbf{c}}(i)$$, and polynomial parameters $$b_{\varphi }$$ that define the axial vertebral rotation $$\varphi (i)$$ in the given 3D spine image.

Moreover, it is assumed that the 3D image is cropped to a volume of interest according to the start and end point of observation along axis $$z$$ of the image-based coordinate system, so that the resulting cropped 3D image contains only axial cross-sections that display the observed anatomy of the spine. Although the presented examples may be therefore labeled as semi-automated, manual determination of the volume of interest does not represent a very demanding or time-consuming task. An advantage of such an assumption is that the parametrization of the spine curve and axial vertebral rotation can be based on axial pixel coordinates $$z$$ at the start and end point of observation, resulting in $$i_{\text{sp}} = 1$$ and $$i_{\text{ep}} = Z$$, respectively, with the corresponding number of samples equal to $$N = i_{\text{ep}} - i_{\text{sp}} + 1 = Z$$, where $$Z$$ is the number of axial cross-sections in the 3D image. Such parametrization is, considering the usual in-plane resolution and slice thickness of CT and MR spine images, in general sufficient for a smooth and continuous description of the spine curve $${\mathbf{c}}(i)$$ and axial vertebral rotation $$\varphi (i)$$.


Automated Determination of the Spine Curve and Axial Vertebral Rotation in CT Images

Vrtovec et al. [83] proposed a method for automated determination of the spine curve and axial vertebral rotation in CT images. If the spine curve $${\mathbf{c}}(i)$$ is represented by a curve that passes through the centers of vertebral bodies, then its determination can be based on the anatomical property that vertebral bodies are locally the largest bone structures of the spine, and on the geometrical property that the center of each vertebral body is represented by the point that is most distant from corresponding edges of the vertebral body. To obtain a quantitative representation of these properties, a distance transform function based on Euclidean metrics is applied twice to image $$I = I({\mathbf{p}})$$, resulting in distance map $$D_{I} = D_{I} ({\mathbf{p}})$$:


$$D_{I} ({\mathbf{p}}) = \left\{ {\begin{array}{*{20}c} { + d_{ < } ({\mathbf{p}},{\tilde{\mathbf{p}}});} & {I({\mathbf{p}}) \ge T,} & {I({\tilde{\mathbf{p}}}) < T,} \\ { - d_{ \ge } ({\mathbf{p}},{\tilde{\mathbf{p}}});} & {I({\mathbf{p}}) < T,} & {I({\tilde{\mathbf{p}}}) \ge T,} \\ \end{array} } \right.$$

(28)
where $$d_{ < } ({\mathbf{p}},{\tilde{\mathbf{p}}})$$ and $$d_{ \ge } ({\mathbf{p}},{\tilde{\mathbf{p}}})$$ are the Euclidean distances between the observed point $${\mathbf{p}} = (x,y,z)$$ and point $${\tilde{\mathbf{p}}} = (\tilde{x},\tilde{y},\tilde{z})$$, which represents the closest point to $${\mathbf{p}}$$ with $$I({\tilde{\mathbf{p}}}) < T$$ and $$I({\tilde{\mathbf{p}}}) \ge T$$, respectively. The image intensity threshold $$T$$, which separates the bone structures from the background, can be in the case of CT images determined from the corresponding Hounsfield values. Each value at point $${\mathbf{p}}$$ in the resulting distance map $$D_{I}$$, which is of the same size as image $$I$$, represents the Euclidean distance from $${\mathbf{p}}$$ to the edges of the bone structures, and this distance is positive when $${\mathbf{p}}$$ is located inside and negative when $${\mathbf{p}}$$ is located outside the bone structures (Fig. 9). As vertebral bodies are locally the largest bone structures of the spine, distance map values are expected to be the highest in geometrical centers of vertebral bodies and smoothly decrease by moving away from the centers. The optimal polynomial parameters $$b_{{\mathbf{c}}}^{ * }$$ that define the spine curve $${\mathbf{c}}(i)$$ (Eq. 25) are finally obtained in an optimization procedure that searches for the combination of parameters that corresponds to the maximal sum of distance map values along the spine curve:

A312884_1_En_8_Fig9_HTML.jpg


Fig. 9
Original image $$I$$ (left) and the corresponding distance transform $$D_{I}$$ with superimposed spine boundaries (right), displayed for a selected a sagittal, b coronal and c axial cross-section of a 3D CT image of a scoliotic spine. In the distance transform $$D_{I}$$, brighter elements represent positive distances, while darker elements represent negative distances



$$b_{{\mathbf{c}}}^{ * } = \mathop {\arg \hbox{max} }\limits_{{b_{{\mathbf{c}}} }} \left( {\sum\limits_{i = 1}^{N} D_{I} ({\mathbf{c}}(i)|b_{{\mathbf{c}}} )} \right).$$

(29)

In order to increase the optimization robustness, two additional mechanisms can be applied. First, the amount of information taken into account during optimization can be increased by considering distance map values within a radius from the spine curve, which can be defined from the quantitative morphometrical vertebral analysis [52–54]. Second, optimization can be designed hierarchically and performed on multiple levels. On the first level, the polynomial functions are initialized with degrees of $$K_{x} = K_{y} = K_{z} = 1$$, i.e. as straight lines. When the optimization reaches the termination criterion, polynomial degrees are increased by $$1$$ and the optimization is restarted on the next level, using the polynomial parameters from the previous level for initialization (optionally, the degree of the polynomial function $$c_{z} (i)$$ can be fixed to $$K_{z} = 1$$, therefore always representing a straight line that describes the longitudinal axis of the spine). The maximal polynomial degree can be determined from the flexion points in normal and scoliotic spinal curvatures, as the number of flexion points of a polynomial function is equal to its degree decreased by one (e.g. a straight line is of the first degree and has zero flexion points). For example, when observing the thoracolumbar section of the spine, three distinctive flexion points exist in a normal spinal curvature, i.e. the maximal thoracic kyphosis, the thoracolumbar junction, and the maximal lumbar lordosis, and therefore polynomial functions of the fourth degree represent an adequate choice. Accordingly, a scoliotic spinal curvature can be adequately described by polynomial functions of the fifth degree.

As the axial vertebral rotation is defined as the rotation of vertebrae about the spine curve, it can be determined automatically in planes that are orthogonal to the spine curve $${\mathbf{c}}(i)$$, defined by polynomial parameters $$b_{{\mathbf{c}}}^{ * }$$ (Eq. 29). In each such $$i$$th plane, a circular domain of radius $$d$$ and centered at point $${\mathbf{c}}(i) = (c_{x} (i),c_{y} (i),c_{z} (i))$$ is defined (Fig. 10a). The circular domain is then bisected by a line that passes through the center of the domain and is inclined for angle $$\varphi$$ against the projection $${\tilde{\mathbf{e}}}_{Iy} (i)$$ of the unit vector $${\hat{\mathbf{e}}}_{Iy} = [0,1,0]_{I}$$ to the plane (Eq. 14). In the obtained halves $$A$$ and $$B$$ of the $$i$$th circular domain, $$s_{A}$$ and $$s_{B}$$ represent image intensities at mirror pixels according to the line of bisection:

A312884_1_En_8_Fig10_HTML.gif


Fig. 10
a A circular domain of radius d and centered at the point on the spine curve $${\mathbf{c}}(i)$$, defined in the plane orthogonal to the spine curve $${\mathbf{c}}(i)$$. b For a bisecting line that is inclined for an angle $$\varphi$$ against the direction of $${\tilde{\mathbf{e}}}_{Iy}$$, image intensities $$s_{A}$$ and $$s_{B}$$ at mirror point pairs are used to evaluate the similarity between the two halves of the circular domain



$$\begin{aligned} s_{A} (i,j,\varphi ) & = I(R_{z} (\varphi ) [ - u,v,c_{z} (i)]) \\ & = I(R_{{{\hat{\mathbf{t}}}(i)}} (\varphi ) [ - x,y,0] + [c_{x} (i),c_{y} (i),c_{z} (i)]), \\ \end{aligned}$$

(30)



$$\begin{aligned} s_{B} (i,j,\varphi ) & = I(R_{z} (\varphi ) [ + u,v,c_{z} (i)]) \\ & = I(R_{{{\hat{\mathbf{t}}}(i)}} (\varphi ) [ + x,y,0] + [c_{x} (i),c_{y} (i),c_{z} (i)]), \\ \end{aligned}$$

(31)
where $$j$$ is the index of the mirror point pair (a total of $$J$$ mirror point pairs exist), consecutively assigned on the basis of each $$(u,v)$$ with $$u > 0$$” src=”/wp-content/uploads/2016/12/A312884_1_En_8_Chapter_IEq237.gif”></SPAN> and <SPAN id=IEq238 class=InlineEquation><IMG alt=, or each $$(x,y)$$ with $$x > 0$$” src=”/wp-content/uploads/2016/12/A312884_1_En_8_Chapter_IEq240.gif”></SPAN> and <SPAN id=IEq241 class=InlineEquation><IMG alt= (Fig. 10b). Matrices $$R_{z}$$ (Eq. 17) and $$R_{{{\hat{\mathbf{t}}}(i)}}$$ (Eq. 19) represent, respectively, the rotation1 about axis $$w$$ of the spine-based coordinate system and the rotation about the unit tangent vector $${\hat{\mathbf{t}}}(i)$$ to the spine curve $${\mathbf{c}}(i)$$ at point $$i$$ in the image-based coordinate system.

Radius $$d$$ is defined so that the anatomy of the whole vertebra (and not only of the vertebral body) is captured within the circular domain, and can be defined from the quantitative morphometrical vertebral analysis [52–54]. From the obtained mirror image intensity pairs, the in-plane similarity between the two mirror halves of the $$i$$th circular domain at inclination $$\varphi$$ can be quantitatively evaluated by the correlation coefficient $$R_{AB} (i,\varphi )$$:


$$R_{AB} (i,\varphi ) = \frac{{\sum\limits_{j = 1}^{J} (s_{A} (i,j,\varphi ) - \bar{s}_{A} (i,\varphi ))(s_{B} (i,j,\varphi ) - \bar{s}_{B} (i,\varphi ))}}{{\sqrt {\sum\limits_{j = 1}^{J} (s_{A} (i,j,\varphi ) - \bar{s}_{A} (i,\varphi ))^{2} \sum\limits_{j = 1}^{J} (s_{B} (i,j,\varphi ) - \bar{s}_{B} (i,\varphi ))^{2} } }},$$

(32)
where $$\bar{s}_{A} (i,\varphi ) = \frac{1}{J}\sum\nolimits_{j = 1}^{J} s_{A} (i,j,\varphi )$$ and $$\bar{s}_{B} (i,\varphi ) = \frac{1}{J}\sum\nolimits_{j = 1}^{J} s_{B} (i,j,\varphi )$$ are the mean image intensities in parts $$A$$ and $$B$$, respectively, of the circular domain. It is important to note that features other than image intensities can be extracted at mirror points (e.g. image intensity gradients), and similarity measures other than the correlation coefficient can be computed (e.g. mutual information) for the corresponding circular domain halves. Nevertheless, the domain must be circular so that irrespectively of inclination $$\varphi$$, the same points are taken into account for the computation of the in-plane similarity. The optimal polynomial parameters $$b_{\varphi }^{ * }$$ that define the axial vertebral rotation $$\varphi (i)$$ (Eq. 26) are finally obtained in an optimization procedure that searches for the combination of parameters that corresponds to the maximal sum of correlation coefficients along the spine curve:


$$b_{\varphi }^{ * } = \mathop {\arg \hbox{max} }\limits_{{b_{\varphi } }} \left( {\sum\limits_{i = 1}^{N} R_{AB} (i,\varphi )|b_{\varphi } } \right).$$

(33)

Similarly as in the case of the spine curve, optimization can be designed hierarchically and performed on multiple levels. However, in the case of the axial vertebral rotation, on the first level the optimization starts with a zero degree polynomial, i.e. $$K_{\varphi } = 0$$, and with polynomial parameter $$b_{\varphi ,0} = 0$$, meaning that the initial axial vertebral rotation is constant and equal to zero along the whole length of the spine. When the optimization reaches the termination criterion, the polynomial degree is increased by $$1$$ and the optimization is restarted on the next level, using the polynomial parameters from the previous level for initialization. The maximal polynomial degree of $$5$$ is adequate for modeling the axial vertebral rotation in both normal and scoliotic spines.

The method was evaluated [83, 85] on $$30$$ normal and one scoliotic CT image of the thoracolumbar spine, and the reported mean difference between the obtained spine curve in 3D and manually defined ground truth points was $$2.1 \pm 1.4$$ mm. The performance of the determination of the axial vertebral rotation was, using the sum of absolute differences instead of the correlation coefficient (Eq. 32), evaluated [18] on $$68$$ vertebrae extracted from CT images, and the reported mean difference against reference values was around $$- 0.6^{ \circ }$$ with the corresponding 95% confidence interval of around $$4.8^{ \circ }$$.


Automated Determination of the Spine Curve and Axial Vertebral Rotation in MR Images

Vrtovec et al. [86] proposed a method for automated determination of the spine curve and axial vertebral rotation in MR images. The method is based on anatomical properties that vertebral bodies and intervertebral discs are nearly circular in shape and that vertebrae are nearly symmetrical over the lines that pass through the centers of vertebral bodies. In each axial cross-section $$z = z_{i}$$ of the MR image (Fig. 11a), an arbitrary in-plane line $$y_{i} = y_{i} (x)$$ that divides the observed cross-section into two image parts, i.e. part $$A$$ and part $$B$$, can be defined as:

A312884_1_En_8_Fig11_HTML.gif


Fig. 11
a The search for the line of in-plane symmetry is performed in each ith axial orthogonal multi-planar cross-section of the MR spine image. b The in-plane line $$y_{i} (x)$$ is defined by parameters $$\gamma_{i}$$ and $$\lambda_{i}$$. c The operator $$\varGamma$$ is centered in $$(x_{i} ,y_{i} )$$ and consists of M concentric rings; $$m = 1,2, \ldots ,M$$, each with radial width of $$\Delta r$$. d The resulting center of the vertebral body $$(x_{i}^{ * } ,y_{i}^{ * } )$$ and in-plane rotation $$\gamma_{i}^{ * }$$



$$y_{i} (x) = \left( {x - \lambda_{i} } \right)\tan \left( {\frac{\pi }{2} - \gamma_{i} } \right),$$

(34)
where $$\tan \left( {\frac{\pi }{2} - \gamma_{i} } \right)$$ is the slope and $$\lambda_{i}$$ is the intersection of the line with axis $$x$$ of the image-based coordinate system (Fig. 11b). The line that splits the observed cross-section into two symmetrical parts can be obtained by maximizing the similarity between image parts $$A$$ and $$B$$:


$$\{ \gamma_{i}^{ * } ,\lambda_{i}^{ * } \} = \mathop {\arg \hbox{max} }\limits_{{\{ \gamma_{i} ,\lambda_{i} \} }} (S(I({\mathbf{p}}_{A} ,z_{i} ),I({\mathbf{p}}_{B} ,z_{i} ))),$$

(35)
where $${\mathbf{p}}_{A} = (x_{A} ,y_{A} );\forall {\mathbf{p}} = (x,y) > y_{i} (x)$$” src=”/wp-content/uploads/2016/12/A312884_1_En_8_Chapter_IEq281.gif”></SPAN> and <SPAN id=IEq282 class=InlineEquation><IMG alt= represent all points $${\mathbf{p}} = (x,y)$$ that lie in image parts $$A$$ and $$B$$, respectively, and $$S$$ is the similarity measure that can be computed as the standard mutual information between image parts $$A$$ and $$B$$:


$$S(I({\mathbf{p}}_{A} ,z_{i} ),I({\mathbf{p}}_{B} ,z_{i} )) = \sum\limits_{{{\mathbf{p}}_{A} \in A}} \sum\limits_{{{\mathbf{p}}_{B} \in B}} {p}_{Q} ({\mathbf{p}}_{A} ,{\mathbf{p}}_{B} )\log \left( {\frac{{{p}_{Q} ({\mathbf{p}}_{A} ,{\mathbf{p}}_{B} )}}{{{p}_{Q} ({\mathbf{p}}_{A} ) {p}_{Q} ({\mathbf{p}}_{B} )}}} \right),$$

(36)
where $$p_{Q} ({\mathbf{p}}_{A} )$$ and $$p_{Q} ({\mathbf{p}}_{B} )$$ are the probability distributions of image intensities in image parts $$A$$ and $$B$$, respectively, $$p_{Q} ({\mathbf{p}}_{A} ,{\mathbf{p}}_{B} )$$ is the joint probability distribution of image intensities, and $$Q$$ is the number of bins used for probability estimation. The resulting parameters $$\gamma_{i}^{ * }$$ and $$\lambda_{i}^{ * }$$ (Eq. 34) define the in-plane line of symmetry, which passes through the center of the vertebral body in the observed $$i$$th axial cross-section.

To determine the exact location of the center of the vertebral body, the shape properties of the vertebral anatomy are combined with the appearance of the vertebral body in MR images. When observed in axial cross-sections, the shape of the vertebral body is relatively circular. For a circular structure displayed in a MR image, a certain amount of variation in image intensities is always present along any radial direction, however, in the tangential direction, the variation in image intensities is relatively small. To obtain a quantitative estimation of these properties, the entropy-based operator $$\varGamma$$ is introduced that is centered at an arbitrary point $${\mathbf{p}} = (x,y)$$ and consists of $$M$$ concentric rings of radii $$\{ r_{m} ; m = 1,2, \ldots ,M\}$$ with $$\forall m {:} \, r_{m} < r_{m + 1}$$ and radial width of each ring equal to $$\Delta r = r_{m + 1} - r_{m}$$ (Fig. 11c):
Tags:
Dec 23, 2016 | Posted by in NEUROLOGICAL IMAGING | Comments Off on Determination of the Spine-Based Coordinate System for an Efficient Cross-Sectional Visualization of 3D Spine Images

Full access? Get Clinical Tree
Get Clinical Tree app for offline access