The potential of MRI, for tracking structural and functional spinal cord changes in vivo in neurological disorders, is currently not thoroughly exploited, both for research and clinically. This is because of many unresolved technical challenges, which involve both image acquisition and processing ¹. With regard to the processing phase, difficulties arise primarily from the small cross-sectional dimensions of the cord, combined with low gray-to-white matter contrast and high morphological variability across subjects. Consequently, as opposed to the large availability of brain image processing software, spinal cord imaging suffers from the lack of similar computational tools. Extracting the cord and segmenting it into gray and white matter represent pivotal steps for most spinal cord MRI studies. Such tasks are typically performed manually, by trained experts, thus making the process time consuming, labor-intensive and fairly subjective. For these reasons, the development of spinal cord segmentation algorithms constitutes a primary research goal in the field of medical image analysis. Nevertheless, the methods that have been proposed so far are mostly semi-automated and require user input, often in the form of a partial labelling or delineation ². Instead, in this work, we present a fully automated and unsupervised framework that can be effectively applied to create probabilistic tissue templates of the spinal cord and segment individual images, with negligible user effort.

Our method exploits probabilistic clustering techniques to learn tissue specific intensity distributions, given MRI scans of a sufficiently large population. Morphological variability is captured simultaneously across subjects in the form of average-shaped tissue templates. Such a representation of the data provides a powerful framework to infer tissue labels without relying on the availability of preexisting atlases. We model the observed image intensities as random variables drawn from Gaussian mixture distributions, with the incorporation of unknown tissue priors that are warped non-linearly to match the individual scans ³. As a result, for a dataset consisting of M individual scans, each one comprising N voxels, the likelihood of $$$\mathbf{X}$$$ can be expressed as

$$p(\mathbf{X}| \mathbf{μ,Σ,π,α})=∏_i^M∏_j^N∑_k^K π_{jk} (α_{ij} )\;\mathcal{N}(x_{ij} |μ_{ik},Σ_{ik})\;,$$

with $$$\{\mathbf{μ},\mathbf{Σ}\}$$$ denoting the Gaussian means and covariances for the K tissue classes, $$$\mathbf{π}$$$ representing the tissue priors and $$$\mathbf{α}$$$ the deformation parameters to map between the individual and the reference anatomical space. Additionally we introduce priors over $$$\mathbf{α}$$$ to assure biophysical plausibility of the deformations and topology preservation. Thus we obtain

$$p(\mathbf{X,α}|\mathbf{μ,Σ,π})=p(\mathbf{α})∏_i^M∏_j^N∑_k^K π_{jk} (α_{ij} )\;\mathcal{N}(x_{ij} |μ_{ik},Σ_{ik})\;.$$

We estimate $$$\{\mathbf{μ},\mathbf{Σ},\mathbf{π}\}$$$ making use of the Expectation-Maximization algorithm for Maximum Likelihood estimation, while we adopt a Gauss-Newton optimization scheme to obtain Maximum a Posteriori estimates of $$$\mathbf{α}$$$. The output of our algorithm consist of: a set of point estimates of the model parameters, segmentations of the individual scans and average shaped, population specific tissue probability maps.

We tested our method on a dataset consisting of 22 scans obtained in healthy control subjects acquired with a 3D high-resolution optimized proton density weighted multi-echo sequence (MEDIC). High-resolution axial volumes of the cervical cord at C2/C3 level were obtained with a resolution of 0.25×0.25×2.50 mm³. The following parameters were used: field of view (FOV) of 162×192 mm², repetition time (TR) of 44 ms, echo time (TE) of 19 ms and flip angle α=11. All subjects are healthy adults (aged between 24 and 56 years). We applied our algorithm by implementing a multiresolution (coarse to fine) strategy consisting of 5 stages with axial resolutions of 1×1mm², 0.8×0.8 mm², 0.5×0.5 mm², 0.4×0.4 mm² and 0.25×0.25 mm² respectively. We set the number of classes equal to 6. Before performing the last two optimization stages we manually restricted the field of view to 37.5×35 mm² so as to focus on the cord area only. The resulting white and gray matter templates are illustrated in figure 1 (b-c) together with an average shaped intensity template (a). Three individual segmentations are reported in figure 2. We used our automatically generated segmentations to calculate the mean spinal cord area, white matter area and gray matter area for all subjects of the population. Results are reported in figure 3, where they are compared to equivalent measurements obtained by expert manual delineation.We presented a general framework to estimate tissue templates from MRI datasets and perform intensity based image segmentation. We demonstrated that the method can be successfully applied to delineate internal structures of the spinal cord and that the segmentations can be used to obtain morphological measurements that are consistent with the ones produced by manual annotation. Therefore our method could represent a fast and reliable alternative to manual spinal cord segmentation.