Semi-automatic Method for the Quantification of Spinal Cord Atrophy

); 192 slices in sagittal orientation were acquired with an in-slice resolution of $0.98\,\text {mm}\times 0.98~\text {mm}$ and a slice thickness of 1 mm. The phantom was scanned in 11 different -positions relative to the magnetic field center ( to 50 mm in increments of 10 mm). The manufacturer’s three-dimensional distortion correction routine was applied to the data to see the effects of image distortions induced by gradient non-linearity [7]. Both original and corrected datasets were reconstructed.

To assess the scan–rescan reliability, 12 healthy volunteers (3 female, 9 male, mean age 32.4 y, range 26–44 y) were scanned on a 3T whole-body MR scanner (Verio, Siemens Medical, Germany) with a T1-weighted MPRAGE sequence (TR/TI/TE/ $\alpha = 2.0\,\mathrm {s}\slash {} 1.0\,\mathrm {s}\slash {}3.4\,\text {ms}\slash {}8^{\circ }$ ); 192 slices in sagittal orientation parallel to the interhemispheric fissure were acquired with an isotropic resolution of $1\,\text {mm}^{3}$ . Both original and distortion-corrected datasets were reconstructed.

To show the applicability to clinical data, 12 relapsing-remitting MS patients (8 female, 4 male, mean age 32.2 y, range 21–46 y; mean disease duration 8.2 y, range 1–17 y, median EDSS 3.0, range 1–4) and 12 age-matched controls (6 female, 6 male, mean age 31.6 y, range 22–48 y) were scanned on a 3T whole-body MR scanner (Verio, Siemens Medical, Germany) with a T1-weighted MPRAGE sequence (TR/TI/TE/ $\alpha =1.6~\mathrm {s}\slash {} 0.9\,\mathrm {s}\slash {} 2.7~\mathrm {ms}\slash {} 9^{\circ }$ ); 192 slices in sagittal orientation parallel to the interhemispheric fissure were acquired with an isotropic resolution of $1~{\text {mm}}^{3}$ . Distortion-corrected datasets were reconstructed.

3 Method

The proposed method can be broken down into four distinct steps, which we refer to as presegmentation, segmentation refinement, surface reconstruction, and reformatting. Of these steps, only presegmentation and reformatting need manual intervention while the others run in a completely automated manner. In this way, the user interaction time lies in the order of two to five minutes per scan.

Fig. 1

a Location of the cisterna pontis. b Presegmentation (blue region of interest, red background seeds, green cord seeds, yellow result). c $$+$$

d Refinement: c intensity profile locations, d single intensity profile (solid) with smooth derivative estimate (dashed), both normalized for display

3.1 Presegmentation

The aim of the presegmentation step (see Fig. 1b) is to gain a binary voxel mask that roughly separates the SC section of interest from the background; that is, from the surrounding cerebrospinal fluid (CSF) and all non-cord tissue. While in principle any kind of thresholding technique could be applied, we use graph cuts [3] because of their flexibility and speed. To compensate for intensity differences caused by field inhomogeneities, we apply a bias field correction [11] to the image volumes beforehand. Furthermore, we normalize the image intensities to the $$[0,1]$$

interval. We then build a six-connected graph from the voxels around the region of interest (which the user may sketch in transverse, sagittal, and coronal projections of the image volume).

The t-link weights are calculated based on a naive Bayes classifier via the intensity distributions of a set of foreground and background seed points; that is, a selection of voxels labeled by the user as definitely belonging either to the SC or its surroundings. We model the foreground as a univariate normal distribution and the background as a mixture of four Gaussians. More specifically, we calculate the weights $w_{\mathrm {fg}}(x)$ and $w_{\mathrm {bg}}(x)$ for the t-links that connect voxel $$x$$

to the foreground and background terminal, respectively, as

$w_{\mathrm {fg}}(x)={\left\{ \begin{array}{ll} \infty &{} \text {if }x\in \fancyscript{F}\\ 0 &{} \text {if }x\in \fancyscript{B}\\ 1-w_{\mathrm {bg}}(x) &{} \text {else} \end{array}\right. }\quad \text {and}\quad w_{\mathrm {bg}}(x)={\left\{ \begin{array}{ll} 0 &{} \text {if }x\in \fancyscript{F}\\ \infty &{} \text {if }x\in \fancyscript{B}\\ \frac{p_{\mathrm {bg}}(I(x))}{p_{\mathrm {bg}}(I(x))+p_{\mathrm {fg}}(I(x))} &{} \mathrm {else} \end{array}\right. },$

where $\fancyscript{F}$ and $\fancyscript{B}$ are the sets of foreground and background seed points, $$I(x)$$

is the intensity of voxel $$x$$

, and $p_{\mathrm {fg}}$ and $p_{\mathrm {bg}}$ are the probability density functions that we estimated from the foreground and background seed point intensities.

The n-link weights $w(x_{a},x_{b})$ between neighboring voxels $x_{a}$ and $x_{b}$ are calculated as $w(x_{a},x_{b})=\kappa \exp (-0.5\varsigma ^{-2}(I(x_{a})-I(x_{b}))^{2})$ , where $I(x_{a})$ and $I(x_{b})$ are the respective voxel intensities, $\kappa$ is a weighting factor, and $\varsigma$ determines the spread of the Gaussian-shaped function (with smaller values for $\varsigma$ leading to a faster decrease of $w(x_{a},x_{b})$ for increasing differences $I(x_{a})-I(x_{b})$ ). The presegmentation is concluded by connected-component labeling, assuring that only the region that includes foreground seeds is retained.

3.2 Segmentation Refinement

Let $I{:}\,\varOmega \rightarrow \mathbb {R}$ denote the preprocessed (i.e., bias field corrected) image, and let $\mathbf {x}=(x_{1},x_{2},x_{3})^{\mathsf {T}}\in \varOmega$ denote the voxel indices in the image domain $\varOmega \subset \mathbb {N}^{3}$ . Furthermore, let $M\subset \varOmega$ denote the set of foreground mask voxel indices from the presegmentation step. To reduce noise in $$I$$

, we apply the GradientAnisotropicDiffusionImageFilter of ITK¹ with the conductance parameter fixed to $$3.0$$

, the time step fixed to $$0.05$$

, and the number of iterations fixed to $$20$$

and

for the image used in the first and second pass, respectively. In this way, we yield the denoised images $\hat{I}_{1}(\mathbf {x})$ and $\hat{I}_{2}(\mathbf {x})$ .

In the first pass, for each transversal slice, we determine the mask boundary voxel indices $B_{z}$ as

$\begin{aligned} B_{z}=M_{z}\backslash (M_{z}\ominus S_{4}), \end{aligned}$

(1)

where $\ominus$ denotes the morphological erosion operator, $S_{4}$ is the two-dimensional structuring element representing four-connectivity, and $M_{z}=\{(x_{1},x_{2},x_{3})^{\mathsf {T}}\in M{:}\,x_{3}=z\}$ is the subset of mask voxels for the $$z$$

-th transversal slice.

We then fit a periodic smoothing B-spline [5] $\mathbf {s}_{z}(\tau )$ of degree three through the ordered voxel indices $\mathbf {b}_{i}^{z}\in B_{z}$ . We distribute the spline’s knots $t_{i}\in [0,1]$ ( $i=1,\ldots ,\left| B_{z}\right|$ ) according to

$\begin{aligned} t_{i}=\tfrac{t_{i}^{*}}{t_{\left| B_{z}\right| }^{*}}\quad \text {with}\quad t_{1}^{*}=0\quad \text {and}\quad t_{i}^{*}=t_{i-1}^{*}+\Vert \mathbf {b}_{i}^{z}-\mathbf {b}_{i-1}^{z}\Vert , \end{aligned}$

(2)

and constrain the spline smoothness by a smoothing parameter $$s$$

via

$\begin{aligned} \Vert \mathbf {b}_{i}^{z}-\mathbf {s}_{z}(t_{i})\Vert ^{2}\le s. \end{aligned}$

(3)

The order of the boundary voxels $\mathbf {b}_{i}^{z}$ is determined by calculating an estimate of their centroid as ${\hat{\mathbf {c}}}_{z}=\frac{1}{\left| B_{z}\right| }\sum _{i}\mathbf {b}_{i}^{z}$ and then sorting them according to the angles formed by the $x_{1}$ -axis and the vectors $\mathbf {b}_{i}^{z}-{\hat{\mathbf {c}}}_{z}$ . Once we have $\mathbf {s}_{z}(\tau )$ , we divide it into $n_{1}$ sections of equal arc length, yielding $n_{1}$ new vertices $\mathbf {u}_{j}^{z}\in \mathbf {s}_{z}(\tau )$ ( $j=1,\ldots ,n_{1}$ ) at the section endpoints. For each $\mathbf {u}_{j}^{z}$ , we then extract a one-dimensional intensity profile (see Fig. 1c+d) $P_{j}^{z}(x)$ with $x\in [0\ldots k_{1}-1]$ as

$\begin{aligned} P_{j}^{z}(x)=\hat{I}_{1}\mathopen {}\left( \mathbf {v}_{j}^{z}+\delta (x)\right) \quad \text {with}\quad \delta (x)=d_{1}\cdot \left( x-\tfrac{k_{1}-1}{2}\right) \mathbf {n}_{j}^{z}\quad \text {and}\quad \mathbf {v}_{j}^{z}=\mathbf {u}_{j}^{z}-o\mathbf {n}_{j}^{z}, \end{aligned}$

(4)

where $\mathbf {n}_{j}^{z}$ denotes the unit normal vector pointing inside the spline curve at $\mathbf {u}_{j}^{z}$ , the resampling distance is given via $d_{1}\in \mathbb {R}$ , the number of profile samples via $k_{1}\in \mathbb {\mathbb {N}}$ , and $o\in \mathbb {R}$ is an offset to control the profile centering with respect to $\mathbf {u}_{j}^{z}$ , yielding offset-corrected vertices $\mathbf {v}_{j}^{z}$ . As our approach to calculate $B_{z}$ systematically underestimates the mask boundary by half a voxel, we set $$o=0.5$$

. To get the intensity values at the resampling positions $$x$$

, we use bilinear interpolation in the $$z$$

th transversal slice of $\hat{I}_{1}$ .

Given that the extracted profiles point to the inside of the spline curve and knowing that in T1-weighted images the inside (i.e., the SC tissue) typically appears brighter than its immediate surroundings (i.e., the CSF), we try to refine the $\mathbf {v}_{j}^{z}$ by means of edge detection in the profiles $P_{j}^{z}$ . We thus calculate their derivatives as $P_{j}^{z\prime }(x)=G^{\prime }(\sigma )*P_{j}^{z}(x)$ , where $G^{\prime }(\sigma )$ is the spatial derivative of a Gaussian kernel with standard deviation $\sigma$ and zero mean. We then search for all local maxima $x_{m,j}$ in $P_{j}^{z\prime }$ and calculate the new boundary estimate ${\mathbf {w}}_{j}^{z}$ as

$\begin{aligned} \mathbf {w}_{j}^{z}=\mathbf {v}_{j}^{z}+\delta (\hat{x}_{m,j})\quad \text {with}\quad \hat{x}_{m,j}= {\arg \min \limits _{x_{m,j}}}\Vert \delta (x_{m,j})\Vert . \end{aligned}$

(5)

To be less susceptible to noise, we dismiss all $x_{m,j}$ with $P_{j}^{z\prime }(x_{m,j})<c\cdot P_{j}^{z\prime }(x_{\mathrm {max},j})$ beforehand, where $c\in [0,1]$ serves as a threshold and $x_{\mathrm {max},j}$ is the global maximum position of $P_{j}^{z\prime }$ . If no valid maxima are retained, which may be the case if the global maximum value is negative, we set $\mathbf {w}_{j}^{z}=\mathbf {v}_{j}^{z}$ .

A second pass of boundary estimation follows, similar to the first pass, starting with the $\mathbf {w}_{j}^{z}$ as initial estimate. The only differences are the following: first, to ensure a homogeneous distribution of the boundary estimates, particularly with regard to the surface reconstruction step (Sect. 3.3), the fitted spline is now resampled $n_{2}$ times at equal angular distances, using the spline center as point of reference, yielding redistributed estimates $\mathbf {y}_{j}^{z}$

Only gold members can continue reading. Log In or Register to continue