the Brain Surface from CT Images with Artifacts Using Dictionary Learning for Non-rigid MR-CT Registration

Fig. 1.

Dictionary-based appearance learning and segmentation of post-surgical cortical brain surface in CT images for direct non-rigid registration of with pre-surgical MRI. We use a training set of post-implantation CT image to learn two models of image appearance, one inside the cortical surface and one outside. We then use these models to segment the cortical surface from test CT images. Using the extracted post-op surface, we perform surface-based registration with the pre-op MRI in order to co-visualize post-implantation electrodes with pre-op imaging data.

2 Methods

From a clinical database of

epilepsy patients at our institution, we have a set of images $\mathcal {I}=\left\{ I_{\text {MR},i}^1,I_{\text {MR},i}^2,I_{\text {CT},i}^2 |i=1,\ldots ,N \right\}$ , where $I_{m,i}^t$ denotes pre-op images acquired at time

and post-op images acquired at time

using imaging modality m for patient i. Following the current practice, for each patient we (i) create pre-op and post-op brain surfaces

and

from the MR images by extracting isosurfaces of the brain masks generated by using the Brain Extraction Tool (BET) [14], (ii) rigidly register $I_{\text {CT},i}^2$ to $I_{\text {MR},i}^2$ to produce the transformation $T_{\text {CT}\rightarrow \text {MR},i}$ by maximizing the normalized mutual information (NMI) similarity metric [16], and (iii) non-rigidly register $I_{\text {MR},i}^2$ to $I_{\text {MR},i}^1$ using a free-form deformation (FFD) [13] with 30 mm B-spline control point spacing and maximizing their NMI to produce the transformation $T_{\text {MR}\leadsto \text {MR},i}$ . We make use of this data to train our model of brain surface appearance in post-surgical CT images.

2.1 Oriented Local Image Appearance

Let $I:\varOmega _I\subset \mathbb {R}^3\mapsto \mathbb {R}$ be the 3D image that maps points from the spatial domain $\varOmega _I$ to image intensity values. We define the set of intensities in an orientable local image patch $\varPhi$ centered about $\mathbf{u}\in \varOmega _I$ :

$\begin{aligned} \varPhi (\mathbf{u},\sigma ) = \left\{ I(\mathbf{u} + \sigma \mathbf{R}\mathbf{t}) |\quad \forall \mathbf{t}\in \varTheta \right\} . \end{aligned}$

(1)

Here, $\mathbf{R}\in \mathbb {R}^{3\times 3}$ is a rotation matrix consisting of a set of orthonormal basis vectors and $\varPhi$ is a set of d image intensity values sampled at patch template points $\varTheta =\left\{ \mathbf{t}_i | \mathbf{t}_i\in \mathbb {R}^3, i=1,\ldots ,d \right\}$ , whose physical size is controlled by a scale term $\sigma \in \mathbb {R}$ . Typically, standard image patches are aligned with the image axes such that $\mathbf{R}$ uses the standard basis $\mathbf{R}=\mathbf{I}$ , the identity matrix, and $\varPhi$ consists of an isotropic grid of sample points centered about the patch origin. For example, a $5\times 5\times 5$ isotropic image patch $\varPhi$ consists of $$d=125$$

sample points $\varTheta$ arranged in a grid about the patch origin.

In this work, we use oriented image patches, where the orientation of the patch $\mathbf{R}$ is determined by the data. Since we are interested in building a model of image appearance both inside and outside the brain surface, we orient image patches according to local surface geometry. For each point on the surface of interest $\mathbf{u}\in S$ , we compute the local surface normal $\mathbf{n}\in \mathbb {R}^3$ and the directions of principal curvature $\mathbf{k}_1,\mathbf{k}_2\in \mathbb {R}^3$ [6]. These normalized vectors form an orthonormal basis $\mathbf{R}=\left[ \mathbf{n}|\mathbf{k}_1 |\mathbf{k}_2 \right]$ with which we orient the patch $\varPhi (\mathbf{u})$ in (1). Figure 2 provides an illustrative comparison between standard and oriented patches at corresponding points along a surface $\mathbf{u}\in S$ . By orienting the patches in this manner, the texture patterns in our patches exhibit greater invariance to changes in location along the surface S.

Fig. 2.

Examples of local oriented image patches compared to their corresponding standard local image patches oriented along the image axes. For the surface S, the local surface normal $\mathbf{n}$ and directions of principal curvature $\mathbf{k}_1,\mathbf{k}_2$ define a local orthonormal basis by which we orient the image patches. We show patch examples from both inside and outside the cortical surface.

2.2 Training the Cortical Surface Appearance Model

To learn our model of brain boundary appearance, we first map the segmented post-op brain surface to the post-op CT image space, $S_i^{2'} = T_{\text {CT}\rightarrow \text {MR},i}^{-1}\circ S_i^2$ , where $\circ$ is the transformation operator. We then create a sparse representation model $\mathbf{D}_\text {in}$ (dictionary) of the intensities inside of the brain and a model $\mathbf{D}_\text {out}$ of the region just outside the brain [9]. These dictionaries will capture the varieties of textural appearance found near the brain surface boundary. We create these models by extracting a training set of overlapping local image patches inside, $\varPhi _{\text {in}}(\mathbf{u}-\alpha \mathbf{n})$ , and outside, $\varPhi _{\text {out}}(\mathbf{u}+\alpha \mathbf{n})$ , points on the surface $\mathbf{u}\in S_i^{2'}$ , where $\mathbf{n}$ is the outward facing local surface normal and $0<\alpha \le 3.0$ mm defines a narrow band region. As in Sect. 2.1, we use the local surface geometry at $\mathbf{u}$ to orient the image patches.

For each appearance class $c=\left\{ \text {in},\text {out}\right\}$ , we create appearance vectors $\mathbf{p}_c\in \mathbb {R}^d$ for each image patch by concatenating the patch values from the set $\varPhi _c$ , where d is the sample dimensionality determined by the chosen patch sampling template $\varTheta$ in Sect. 2.1, and normalize $\mathbf{p}_c$ to have unit length. Then, we model the distribution of $\mathbf{p}_c$ ’s from all N training images using an overcomplete dictionary $\mathbf{D}_c\in \mathbb {R}^{d\times n}$ such that $\mathbf{p}_c \approx \mathbf{D}_c\varvec{\gamma }$ [1]. Here, n is the number of dictionary atoms and $\varvec{\gamma }$ is the sparse dictionary weighting coefficients. To reconstruct a given appearance sample $\mathbf{p}$ (normalized to have unit length) from the dictionary with a given target sparsity constraint $\varGamma _0$ , we solve the sparse coding problem

$\begin{aligned} \mathop {\min }\limits _\gamma \Vert \mathbf{p}-\mathbf{D}\varvec{\gamma } \Vert _2^2 \text { s.t. } \Vert \varvec{\gamma }\Vert _0\le \varGamma _0, \end{aligned}$

(2)

using an orthogonal matching pursuit (OMP) algorithm [12]. Next, we define the residual error

$\begin{aligned} R_c(\mathbf{p}) =\Vert \mathbf{p}-\mathbf{D}_c\varvec{\gamma }\Vert _2 \end{aligned}$

(3)

for both the inside and outside region classes. Using normalized appearance vectors, $0\le R_c(\mathbf{p})\le 1$ , where values of 0 correspond to perfect signal reconstruction (strong membership to class c) and values of 1 indicate $\mathbf{p}$ could not be reconstructed by $\mathbf{D}_c$ (poor membership to class c). Intuitively, $R_\text {in}(\mathbf{p}) < R_\text {out}(\mathbf{p})$ for points inside the cortical surface boundary. In contrast to Huang et al. [9]’s multi-scale modeling where multi-scale data is concatenated into appearance vectors, we perform this sampling and training procedure individually for each level $k=1,\ldots ,K$ in a multi-resolution Gaussian image pyramid, which means that we train K dictionary pairs $\left\{ \mathbf{D}_\text {in}^k,\mathbf{D}_\text {out}^k \right\}$ .

2.3 Segmenting the Cortical Surface in Post-implantation CT Images

Given a new pair of pre-implantation MRI and post-implantation CT test images not in the training set $I_\text {MR}^1,I_\text {CT}^2\notin \mathcal {I}$ , we perform an initial brain surface segmentation estimate $\hat{S}$ using an intensity-based rigid registration of $I_\text {MR}^1$ and $I_\text {CT}^2$ and transforming the segmented MR surface $$S^1$$

to post-op imaging space, i.e. $\hat{S} = T_{\text {MR}\rightarrow \text {CT}}\circ S^1$ . For a point on the estimated surface $\mathbf{u}\in \hat{S}$ , we extract the oriented local image patch $\varPhi (\mathbf{u})$ to create the appearance vector $\mathbf{p}$ , and then compute the difference of the appearance model residuals from (3) such that

$\begin{aligned} D(\mathbf{p}) = R_\text {out}(\mathbf{p}) - R_\text {in}(\mathbf{p}). \end{aligned}$

(4)

Intuitively, if $\mathbf{u}$ lies within the true boundary of the cortical surface in the CT image then $$$D(\mathbf{p})>0$$” src=”/wp-content/uploads/2016/09/A339424_1_En_52_Chapter_IEq75.gif”></SPAN>, and if <SPAN id=IEq76 class=InlineEquation><IMG alt=$$\mathbf{u}$$ src=$ is outside the true boundary then $D(\mathbf{p})<0$ . Thus, the cortical surface boundary is located at the point $\mathbf{u}$ where $D(\mathbf{p})=0$ . We therefore seek the surface that minimizes the objective function

$\begin{aligned} \hat{S} = \mathop {\min }\limits _S \int _S \Vert D\Vert _2 dS. \end{aligned}$