Robust Cochlear Modiolar Axis Detection in CT

TL;DR

This paper introduces a novel, robust algorithm for detecting the cochlear modiolar axis in CT images, improving accuracy over manual methods by using a kinematic shape representation and advanced surface fitting techniques.

Contribution

The paper presents the first use of a kinematic shape representation in medical image analysis for cochlear axis detection, enhancing robustness and reliability.

Findings

Algorithm reduces alignment error compared to manual detection.

Performs well even with partial surface data.

Validated on synthetic and clinical CT datasets.

Abstract

The cochlea, the auditory part of the inner ear, is a spiral-shaped organ with large morphological variability. An individualized assessment of its shape is essential for clinical applications related to tonotopy and cochlear implantation. To unambiguously reference morphological parameters, reliable recognition of the cochlear modiolar axis in computed tomography (CT) images is required. The conventional method introduces measurement uncertainties, as it is based on manually selected and difficult to identify landmarks. Herein, we present an algorithm for robust modiolar axis detection in clinical CT images. We define the modiolar axis as the rotation component of the kinematic spiral motion inherent in the cochlear shape. For surface fitting, we use a compact shape representation in a 7-dimensional kinematic parameter space based on extended Pl\"ucker coordinates. It is the first time such a kinematic representation is used for shape analysis in medical images. Robust surface fitting is achieved with an adapted approximate maximum likelihood method assuming a Student-t distribution, enabling axis detection even in partially available surface data. We verify the algorithm performance on a synthetic data set with cochlear surface subsets. In addition, we perform an experimental study with four experts in 23 human cochlea CT data sets to compare the automated detection with the manually found axes. Axes found from co-registered high resolution micro-CT scans are used for reference. Our experiments show that the algorithm reduces the alignment error providing more reliable modiolar axis detection for clinical and research applications.

Equations

$$d_i(\mathbf{m})=\frac{\mathbf{v}({\mathbf{p}}_i)\cdot {\mathbf{n}}_i}{\sqrt{\Vert \mathbf{v}({\mathbf{p}}_i){\Vert }^2+w_{\mathbf{p}}\Vert {\nabla }_{\mathbf{p}}(\mathbf{v}({\mathbf{p}}_i)\cdot {\mathbf{n}}_i){\Vert }^2}}$$

$$\log p(\mathcal{D}|\mathbf{m})=\log \prod _{i=1}^{n}p(d_i|\mathbf{m})=-\frac{n}{2}\log 2\pi \Sigma -\frac{1}{2}\sum _{i=1}^n\frac{d_i^2}{\Sigma }=-\frac{n}{2}\log 2\pi \Sigma +\mathcal{L}(\mathbf{m}).$$

$${\mathbf{M}}_i=\mathbf{f}({\mathbf{p}}_i,{\mathbf{n}}_i)\mathbf{f}({\mathbf{p}}_i,{\mathbf{n}}_i{)}^T,{\mathbf{N}}_i=\left[\begin{array}{ccc}{\mathbf{A}}_{{\mathbf{p}}_i}^T{\mathbf{A}}_{{\mathbf{p}}_i}+w_p{\mathbf{A}}_{{\mathbf{n}}_i}^T{\mathbf{A}}_{{\mathbf{n}}_i}& -{\mathbf{A}}_{{\mathbf{p}}_i}& 0\\ -{\mathbf{A}}_{{\mathbf{p}}_i}^T& \mathbf{I}& {\mathbf{p}}_i\\ 0& {\mathbf{p}}_i^T& {\mathbf{p}}_i\cdot {\mathbf{p}}_i+w_p\end{array}\right].$$

$${\mathbf{B}}_m\mathbf{m}={\mathbf{C}}_m\mathbf{m}$$

$$-\psi\Big{(}\frac{\nu}{2}\Big{)}+\log\Big{(}\frac{\nu}{2}\Big{)}+1+\psi\Big{(}\frac{\nu+1}{2}\Big{)}-\log\Big{(}\frac{\nu+1}{2}\Big{)}+\frac{1}{n}\sum_{i=1}^{n}(\log z_{i}-z_{i})=0,$$

$$\Delta\theta=\arcsin{\frac{\|\mathbf{r}_{\mathrm{ref}}\times\mathbf{r}\|}{\|\mathbf{r}_{\mathrm{ref}}\|\|\mathbf{r}\|}},\qquad\Delta d=\Bigl{|}(\mathbf{p}_{0}-\mathbf{p}_{\mathrm{ref}})-\frac{\mathbf{r}\cdot(\mathbf{p}_{0}-\mathbf{p}_{\mathrm{ref}})}{\|\mathbf{r}\|}\Bigr{|}$$

Full text

1 Introduction

The cochlea is a spiral structure in the inner ear that transduces acoustic waves into electrical nerve impulses to enable hearing. The morphology of the human cochlea is complex and highly variable. Therefore, an unambiguous description of the morphological parameters for both modeling and clinical applications is required. Of great importance is the modiolar axis, the central axis of the spiral shape, as it is used to define the z-axis of cylindrical cochlear coordinate reference systems [1]. Due to the tonotopic organization of the cochlea, the modiolar axis connects anatomical features with physiological parameters, mapping spatial positions along the spiral with the perceived characteristic frequencies. In cochlear implantation, in which an electrode array is inserted into the cochlea to restore hearing, radiological parameters referenced by the modiolar axis are used for preoperative planning (selection of suitable implant lengths), for postoperative evaluation (array insertion depth assessment) and for audio processor programming (patient-specific tonotopic stimulation maps). They are also used to investigate the effects of tonotopic mismatch between different stimulation channels on speech rehabilitation in bilateral cochlear implant users.

A common definition of the modiolar axis is based on anatomical landmarks that can only be identified imprecisely in computed-tomography (CT) images: the helicotrema and the center of the modiolus in the basal turn of the cochlea [2]. This leads to misalignment and inter-observer variability, even when using multi-planar reconstructions. As a consequence, outcome measures that are referenced by modiolus-based coordinate systems are distorted. Furthermore, misclassification of cochlear morphology can be caused by inaccurate modiolar axis estimation [3]. Previous detection algorithms use center-line based methods, however either requiring fully segmented image data of the cochlea [3] or knowledge of additional extrinsic parameters [4].

Herein, we present a novel approach for modiolar axis detection suitable for clinical resolution CT images. The spiral shape of the cochlea is modelled as a kinematic surface to mimic its natural growth. Kinematic surfaces are defined as the location of surface points that are tangent to a parameterized stationary velocity field. Then, the modiolar axis is determined as the rotation component of the intrinsic spiral motion. Our contribution lies in the first application of a compact seven-dimensional kinematic surface representation for medical image analysis. Furthermore, we extend the method by a robust maximum likelihood scheme based on a Student-t distribution. The algorithm output is verified using a synthetic data set. Finally, we perform an experimental validation study to compare the modiolar axis detection results between the conventional landmark-based method and our algorithm under consideration of µCT reference data.

2 Methods

2.1 Kinematic Modiolar Axis Detection

To find the modiolar axis, we aim to determine the rotation component of the intrinsic kinematic spiral motion forming the cochlea. We base our work on a line element geometry approach for kinematic surface recognition [5]. A kinematic surface is a surface consisting of oriented points (with position $\mathbf{p}_{i}$ and unit surface normals $\mathbf{n}_{i}$) that is tangent to a parametric velocity field $\mathbf{v}(\mathbf{p})$, i.e. $\mathbf{v}(\mathbf{p})\cdot\mathbf{n}=0$. Kinematic surfaces somewhat extend the notion of implicit surfaces $\mathrm{S}(\mathbf{p})=0$ defined on points.

For the cochlea, we consider a spiral velocity field $\mathbf{v}(\mathbf{p})=\mathbf{r}\times\mathbf{p}+\mathbf{c}+\gamma\mathbf{p}$, consisting of a rotation $\mathbf{r}$, translation $\mathbf{c}$ and scale factor $\gamma$ (see Fig. 1). By choosing $\mathbf{f}(\mathbf{p},\mathbf{n})=\{\mathbf{p}\times\mathbf{n},\mathbf{n},\mathbf{p}\cdot\mathbf{n}\}$ the surface is transformed into a seven-dimensional parameter space based on extended Plücker coordinates. Then, the problem of spiral shape recognition reduces to fitting a linear subspace to $\mathbf{f}(\mathbf{p},\mathbf{n})$ [5]. This is achieved by searching for the parameters $\mathbf{m}=\{\mathbf{r},\mathbf{c},\gamma\}$ such that the distance $d_{i}(\mathbf{m})$ between each point and the surface tangent to the velocity field $\mathbf{v}(\mathbf{p}_{i})$ is minimized. By using a first order approximation of the distance (approximate maximum likelihood method [6, 7]) we can write:

[TABLE]

where $w_{p}$ is a scalar regularizing the denominator. For a spiral velocity field, we find $\nabla_{\mathbf{p}}(\mathbf{v}(\mathbf{p}_{i})\cdot\mathbf{n}_{i})=\mathbf{n}_{i}\times\mathbf{r}+\gamma\mathbf{n}_{i}=[\mathbf{A}_{\mathbf{r}}+\gamma\mathbf{I}]\mathbf{n}_{i}$ where $\mathbf{A_{r}}$ is the skew-symmetric matrix associated with vector $\mathbf{r}$, and $\mathbf{I}$ is the identity matrix. Then we have $\|\nabla_{\mathbf{p}}(\mathbf{v}(\mathbf{p}_{i})\cdot\mathbf{n}_{i})\|^{2}=\|\mathbf{n}_{i}\times\mathbf{r}\|^{2}+\gamma^{2}$ and $\|\mathbf{v}(\mathbf{p}_{i})\|^{2}=\|\mathbf{r}\times\mathbf{p}_{i}\|^{2}+\gamma^{2}\|\mathbf{p}_{i}\|^{2}+2\gamma(\mathbf{p}_{i}\cdot\mathbf{c})+2[\mathbf{r},\mathbf{p}_{i},\mathbf{c}]+\|\mathbf{c}\|^{2}$. In a first approach, we assume a Gaussian distribution of the distance error with variance $\Sigma$, i.e., $p(d_{i})={\cal N}(d|0,\Sigma)$. The objective now is to determine $\mathbf{m}$ such that the log-likelihood $p(\mathcal{D}|\mathbf{m})$ is maximized:

[TABLE]

We can write $d_{i}^{2}=\frac{(\mathbf{v}(\mathbf{p}_{i})\cdot\mathbf{n}_{i})^{2}}{\|\mathbf{v}(\mathbf{p}_{i})\|^{2}+w_{p}\|\nabla_{p}(\mathbf{v}(\mathbf{p}_{i})\cdot\mathbf{n}_{i})\|^{2}}=\frac{\mathbf{m}^{T}\mathbf{M}_{i}\mathbf{m}}{\mathbf{m}^{T}\mathbf{N}_{i}\mathbf{m}}$ with the $7\times 7$ matrices $\mathbf{M}_{i}$, $\mathbf{N}_{i}$ defined as:

[TABLE]

Maximizing the log likelihood $\mathcal{L}(\mathbf{m})$ is equivalent to minimizing $\sum_{i=1}^{n}\frac{\mathbf{m}^{T}\mathbf{M}_{i}\mathbf{m}}{\mathbf{m}^{T}\mathbf{N}_{i}\mathbf{m}}$ which leads to solving the generalized eigenvalue problem [7]:

[TABLE]

with $\mathbf{B}_{m}=\sum_{i=1}^{n}\frac{\mathbf{M}_{i}}{\mathbf{m}^{T}\mathbf{N}_{i}\mathbf{m}}$ and $\mathbf{C}_{m}=\sum_{i=1}^{n}\frac{\mathbf{m}^{T}\mathbf{M}_{i}\mathbf{m}}{(\mathbf{m}^{T}\mathbf{N}_{i}\mathbf{m})^{2}}\cdot\mathbf{N}_{i}$.

The non-linear problem can be tackled by iteratively computing the matrices $\mathbf{B}_{m}$ and $\mathbf{C}_{m}$ for a given estimation of $\mathbf{m}$ and then estimating $\mathbf{m}$ as the eigenvector associated with the smallest eigenvalue (closest to zero). We then define the direction of the modiolar axis as the rotation component $\mathbf{r}$ of the estimated parameters. We further need to compute the zero velocity center of the spiral motion $\mathbf{p}_{0}=\frac{1}{\gamma(\mathbf{r}^{2}+\gamma^{2})}(\gamma\mathbf{r}\times\mathbf{c}-\gamma^{2}\mathbf{c}-(\mathbf{r}\cdot\mathbf{c})\mathbf{r})$ to define the position of the modiolar axis.

2.2 Robust Detection

Like any least-squares fitting method, the above mentioned approach is sensitive to outliers. To increase robustness, we propose to replace the Gaussian likelihood with a Student-t distribution, which is a Gaussian Scale Mixture. More precisely, we assume $p(d_{i})={\mathrm{St}}(d_{i}|0,\Sigma,\nu)=\int_{z_{i}}\mathcal{N}(d_{i}|0,\Sigma/z_{i})\,\mathrm{Ga}(z_{i}|\nu/2,\nu/2)$ where $z_{i}$ is the variance scale variable which has a prior given by the Gamma distribution parameterized by the degrees of freedom $\nu$. When $\nu\Longrightarrow+\infty$, then the Student-t is equivalent to the Gaussian distribution. Therefore, the $\nu$ variable is inversely proportional to the number of outliers. The estimation of the kinematic surface is now extended with an Expectation-Maximization scheme, where $z_{i}$ is the latent variable [8]. This is equivalent to iteratively estimating $\nu$, $z_{i}$ and $\Sigma$ with the following steps:

• •

E-step: Estimate $z_{i}$ for each data point as $z_{i}=(\nu+1)/(\nu+d_{i}^{2}/\Sigma)$. When $\nu$ is very large then $z_{i}$ is close to 1, irrespective to the Mahalanobis distance $d_{i}^{2}/\Sigma$. When $\nu$ is less large then $z_{i}$ is close to zero for outliers (since the Mahalanobis becomes large) and close to 1 for inliers.

• •

M-$\Sigma$ step: Estimate the variance as $\Sigma=\frac{1}{n}\sum_{i=1}^{n}z_{i}d_{i}^{2}$

• •

M-$\nu$ step: Estimate $\nu$ as the solution of the non-linear problem

[TABLE]

where $\psi(x)$ denotes the digamma function.

As an output, we get the estimation of confidence $z_{i}$ in each data point. To obtain a robust estimation of the parameter vector $\mathbf{m}$, we proceed as before (2) but with $\mathbf{B}_{m}=\sum_{i=1}^{n}z_{i}\frac{\mathbf{M}_{i}}{\mathbf{m}^{T}\mathbf{N}_{i}\mathbf{m}}$ and $\mathbf{C}_{m}=\sum_{i=1}^{n}z_{i}\frac{\mathbf{m}^{T}\mathbf{M}_{i}\mathbf{m}}{(\mathbf{m}^{T}\mathbf{N}_{i}\mathbf{m})^{2}}\mathbf{N}_{i}$.

2.3 Implementation

The detection algorithm was implemented in Matlab. It is initialized by specifying the landmarks $\mathbf{l}_{1}$ and $\mathbf{l}_{2}$ (Fig. 2 left) in CT to preselect a spherical segment volume approximately covering the middle and apical turns of the cochlea (radius $r_{s}=\|\mathbf{r}_{s}\|=\|\mathbf{l}_{1}-\mathbf{l}_{2}\|$ and center $\mathbf{c}_{s}=(\mathbf{l}_{1}+\mathbf{l}_{2})/2$, cropped by the planes perpendicular to $\mathbf{r}_{s}$ in $\mathbf{l}_{1}$ and $\mathbf{l}_{2}$). This volume is chosen to minimize interference of proximal structures that may appear connected to the cochlea in clinical CT (i.e., the tympanic cavity at the round window, the internal auditory canal, and the facial nerve). The data is labelled through intensity thresholding (isovalue at 1000 HU), smoothed, and isosurfaces are generated by a marching cubes routine. The largest connected surface with the center of gravity closest to $\mathbf{c}_{s}$ is extracted. The point cloud is scaled and centered, surface normals are computed and the parameter space $\mathbf{f}(\mathbf{p})$ is obtained. Kinematic surface fitting is performed by 5 iterations with $w_{p}=0.001$ (Fig. 2 right).

2.4 Verification in Synthetic Data

We used a polynomial cochlea model with known modiolar axis to generate a synthetic data set for algorithm verification [9]. To mimic the facial nerve located close to the cochlea and often causing segmentation artifacts, we added a tubular structure (1 mm in diameter) with constrained random alignment. In addition, point positions were perturbed with Gaussian noise (0.15 mm standard deviation). We tested the robustness of the algorithm in varying levels of surface coverage (analogous to algorithm implementation, however with different radii covering 5% to 100% of total points). For each level, 500 random cochleae were generated. For comparison, the non-robust (Gaussian) version of the detection algorithm as well as simple detection using principal component analysis (PCA) of the point cloud was applied, where the modiolar axis was selected as the component with least variance. We assessed the alignment between the estimated axis $\mathbf{r}$ and reference axis $\mathbf{r}_{\mathrm{ref}}$ with the angular error $\Delta\theta$ and the distance error $\Delta d$ defined as the closest absolute distance between the estimated axis $\mathbf{r}$ and the center point $\mathbf{p}_{\mathrm{ref}}$ on the reference axis $\mathbf{r}_{\mathrm{ref}}$:

[TABLE]

2.5 Experimental Validation

We validated the algorithm with a data set of 23 human temporal bone specimens consisting of clinical CT (voxel size: 156 $\times$ 156 $\times$ 200 µm3) and co-registered µCT (voxel size: 603 µm3) scans. Four experts manually identified 2 landmarks (see Fig. 2) for each sample in multi-planar CT reconstructions to specify the modiolar axis. The same landmarks were also used to initialize the robust detection algorithm. Again, we applied the simple PCA-based and the non-robust axis estimations for comparison. As reference, the modiolar axis and its center were determined in each specimen using high-resolution surface models from the segmented µCT data. The alignment differences were assessed with the equations shown in (3). Differences in alignment errors were estimated using (separate) linear mixed-effects models (R environment with lme4 package) [10], with a fixed effect for the detection method (categorical variable). We included random intercepts for specimens and experts to account for paired measurements. Before analysis, the data was log transformed.

3 Results

Figure 3 illustrates the verification results. With increasing surface data available, the robust algorithm converges to an angular error of 2.5 degrees and an distance error below 0.1 mm. In contrast, the vertex-based PCA detection shows limited improvement and even yields worse distance errors with increasing cochlea surface coverage, since the basal turn vertices cause a shifting away from the modiolar axis. As expected, the non-robust (Gaussian) version of the algorithm is more sensitive to outliers (facial nerve structure).

Figure 4 summarizes the alignment error after manual landmark-based and automated modiolar axis detection in the CT data of the 23 specimens. The PCA based detection performs worse than manual selection. As measured by the linear mixed effects models, compared with the manual procedure, the robust method reduced the average distance error from 0.32 mm to 0.13 mm (improvement by 0.19 mm, 95% confidence interval [0.17 mm, 0.21 mm]) and the angular error from 9.0° to 2.4° (improvement by 6.6°, 95% confidence interval, [6.1°, 6.9°]). The non-robust version performed worse than the robust version (average distance error 0.04 mm higher and angular error 1.5° higher). The robust procedure further reduced the variability between the observers. Figure 5 visualizes an example.

4 Conclusions

We present a novel and anatomically meaningful approach to model the spiral shape of the cochlea as a structure formed by a kinematic motion (natural growth) and extract the modiolar axis as the rotation component. This formalism enables us to detect the modiolar axis even in subsets of cochlea surfaces obtained from CT images. The approach generalizes the implicit surface representation of oriented points. The parameters $\mathbf{m}=\{\mathbf{r},\mathbf{c},\gamma\}$ provide a compact representation of the cochlea based on intrinsic shape properties and enable novel approaches to cochlear morphology classification. Using appropriate curves, it could be used to generate surfaces for cochlea segmentation. The approach is further applicable for kinematic surface detection of cylindrical, conical, rotational and helical motions [5]. It could be extended by improved surface extraction methods or by directly using image gradients from CT for kinematic parameter space computation. The script runs in $\sim$3 sec on a standard laptop (Intel i7).

4.0.1 Acknowledgments.

Supported by the Swiss National Science Foundation (grant no. P400P2_180822) and the French government (UCAJEDI - ANR-15-IDEX-01).