Automatic localization of cochlear implant electrodes using cone beam computed tomography images | BioMedical Engineering OnLine

Dataset

The dataset used for evaluation consists of postoperative CBCT scans from 150 implanted CI users including 10 different EAs (15 each) (see Table 1). The dataset was created to analyze the accuracy of the localization algorithm if applied to various EAs differing in the number and size of the electrodes or the spacing between them. The CBCT scans were acquired either with a Xoran Minicat (South State Road, Ann Arbor, United States) or with a Xoran xCat (South State Road, Ann Arbor, United States). The resolution of all CBCT scans was 0.3 × 0.3 × 0.3 mm (except one has the resolution 0.125 × 0.125 × 0.125 mm). The selected dataset contained EAs from the four CI manufactures MED-EL (Innsbruck, Austria), Advanced Bionics (Valencia, CA, Unites States), Oticon Medical (Vallauris, France) and Cochlear (Macquarie Park, NSW, Australia). CBCT scanning was performed postoperatively for all subjects, following the implantation. Tables 2 and 3 provide an overview of the CT scanners and the corresponding settings used for each CI user. The seven most basal electrodes of the MED-EL Flex EA consist of double contacts, and the five most apical are single electrode contacts. For this reason, the CBCT scans resulted in lower intensities for the apical electrodes compared to the basal electrodes. The MidScala and SlimJ EAs from Advanced Bionics have 16 active contacts and one non-stimulating marker contact. For these EAs, the distance between the individual active electrodes is constant and the distance between the most basal electrode and the marker is 3 mm. The Neuro ZTI EVO EA from Oticon Medical has also a constant distance between the electrodes. The EAs from Cochlear contain 22 electrodes. The Hybrid-L EA from Cochlear also features a non-stimulating marker contact.

For the entire dataset, a clinical expert determined the electrode locations manually in two different ways. For distantly spaced electrodes, locations were determined labeling the center of the individual electrodes. The location of closely spaced electrodes was determined using a spline that was fitted through the EA. Individual electrode locations were determined using the known geometry of the EA according to Krüger et al. [20]. The manually determined electrode locations were used as ground truth (GT).

The dataset was used to evaluate the improvement that the extension of the cost function we proposed has on the detection rate and localization accuracy.

Evaluation criteria

The localization accuracy \(L_{\text{a}}\) and the detection rate were used to evaluate the automatic electrode localization algorithm. The localization accuracy is based on the averaged localization error defined as the Euclidean distance between a predicted electrode \({\text{el}}_{{\text{pd}},N}\) and the GT electrode position \({\text{el}}_{{\text{GT}},N}\) with the same label or electrode number N (see Eq. 1):

The detection rate was defined as the number of correctly detected electrodes \({\text{el}}_{{\text{detect}}}\) divided by the total number of electrodes \({\text{el}}_i\) in the EA (Eq. 2). An electrode is correctly detected if the Euclidian distance between the location of the automatically predicted electrode \({\text{el}}_{{\text{pd}}}\) and the location of the GT electrode \({\text{el}}_{{\text{GT}}}\) was less than or equal to 0.9 mm (3 voxels):

Graph-based algorithm to locate cochlear implant electrodes

In the present study, automated electrode localization was performed using an implementation of the graph-based procedure according to Hachmann and Nogueira [7] which can be considered own implementation of Noble and Dawant [8]. Figure 5 illustrates the processing steps of the algorithm as a block diagram. It takes as input a region of interest (ROI). First, candidate voxels representing potential electrode locations are selected. The next steps consist of a path finding algorithm for a rough electrode localization and a refinement procedure enabling sub-voxel electrode localization.

Block diagram illustrating the processing steps to determine the cochlear implant electrode position. This process follows the same structure as Noble and Dawant [8] and Hachmann and Nogueira [7]

The input for the localization procedure is a ROI obtained from CBCT data containing the electrode contacts. The ROI was determined based on the GT coordinates. The automatic segmentation of the ROI needs to be implemented in future versions of the algorithm. For this purpose, based on the GT coordinates, a region was defined around the EA with a minimum distance of 3 mm (10 voxels) to the electrodes. Next, a binary mask was generated from the ROI defining a threshold. For this purpose, the cumulative histogram of the intensity values of the voxels contained in the ROI were calculated to define a threshold that depends on the ROI’s dynamic range. According to Noble and Dawant [8], 0.08% of the maximum value of the ROI’s cumulative histogram can be used as threshold α₁. From the contours of the resulting binary mask, centerlines were determined whose voxel coordinates served as candidate points for possible electrode locations [17]. From the candidate points, the path finding algorithm generated all possible paths p consisting of L nodes. The number of nodes (L) corresponds to the number of electrodes in the EA. In contrast to Noble and Dawant [8], we used the first node given by the GT corresponding to the most basal electrode location to reduce computing time. According to Hachmann and Nogueira [7], the most basal GT coordinate serves as a reference for the seed node of the path finding algorithm. This significantly shortens the evaluation time while having minimal impact on the determination of localization accuracy [7]. Paths were built up starting from the seed node as first path node p₁. These paths were evaluated using an initial cost function as shown in Eq. (4):

The initial cost function Cost_initial consists of an intensity-based C_I and an initial shape-based C_S,initial component, where c are the selected candidate points and p are the points already located in the path, which in the case of the initial cost function, was the initial seed point. The intensity-based component is calculated as described below (see Eq. (5)):

In Eq. (5), I_max is the maximum intensity of the ROI, I(c) is the intensity of child node c, i is the length of the path p. α₃ and α₄ are set to 0.1 and 14. The initial shape-based component is calculated according to Eq. (6)):

The next step is to sort out paths in which the two nodes in the path do not have the desired distance to each other. The desired distance is obtained from the inter electrode spacing defined by the manufactured specifications. Further child nodes c are added to the path as path node p_i+1 if located within a certain distance from the path as nodes p_i (Eq. 7):

In Eq. (7), d_i is the distance between electrodes with indices i and i + 1 counted from apex to base. In the equation \(\left\| {{{\varvec{p}}}_{{\varvec{i}}} – {{\varvec{c}}}} \right\|_2\) describes the Euclidean distance between two vector points, which is also called norm 2. This is defined as follows: \(\left\| x \right\|_2 = \sqrt {x_1^2 + x_2^2 + x_x^2 }\), with \(x = p_i – c\). Furthermore, a child node should not already exist in the path. The cost of each possible child node is calculated using the cost function shown in Eq. (8). This cost function consists of an intensity-based \(C_{\text{I}}\) and a shape-based part \(C_{\text{S}}\). Equation (8) shows the calculation of intensity-based costs:

In the intensity-based component, child nodes that correspond to a high intensity in the ROI are preferred. For the last \(L + 1 – \alpha_4\) electrodes, this cost term is reduced because apical electrodes are less bright in the CBCT image. Equation (9) shows the calculation of the cost of the shape-based part C_S:

According to Noble and Dawant [8], the constants in Eq. (9) were set to α₅ = 1.0, α₆ = 5.2 and α₇ = 2.0. The first part of the shape-based cost function C_s was a smoothing term that incurs a high cost if adding a child node would result in a sharp bend in the EA. The second part is a distance term that rejects child nodes that do not have the expected distance \(d_i\) to the last node in path p. The calculation of the cost is performed for all possible child nodes and is added to the already existing costs, from the previous iterations (see Eqs. 11, 12):

where

After this the P = 10 paths with lowest cost are saved. This is repeated until the number of nodes L in the path is reached. The path with the lowest cost is selected as the optimal path.

In the last step, the optimal path is refined so that the position of each node can be estimated on sub-voxel positions. First, a rectangular grid is placed around each point of L,\(\left\{ {{{\varvec{l}}}_{{\varvec{i}}} } \right\}_{i = 1}^{{\varvec{L}}}\) to define node points, which are sampled by \(\left\{ {{\varvec{n}}} \right\}^i = \left\{ {{{\varvec{l}}}_i + \alpha_8 \left[ {x, y, z} \right]} \right\}_{x, y, z \in \left[ { – \alpha_9 ,\alpha_9 } \right]}\). Here, \(\alpha_8\) and \(\alpha_9\) are defined to be 0.12 and 3 mm, respectively. This step aims at refining each estimated position \({{\varvec{l}}}_i\) of every ith electrode with a nearby candidate \(\left\{ {{\varvec{n}}} \right\}^i\). The cost is thereby determined by means of a second cost function (see Eq. (13)):

Here, \(G_\sigma \left( {I\left( {{\varvec{c}}} \right)} \right)\) corresponds to the Gaussian filter response of the image with \(\sigma = 0.3\,{\text{mm}}\). \(\alpha_{10}\) and \(\alpha_{11}\) are set to 50 and 20, respectively. The first term of the second cost function is a scaled blob finding filter. The second term penalizes child nodes c that are not within the expected distance from the last found path electrode p_i. The path finding process is the same as the previous localization step.

New advances and parametrization for automatic localization of cochlear implant electrodes

As reported by Hachmann and Nogueira [7], the Noble and Dawant [8] method was not able to correctly detect the position of all electrodes. First, it was found that double detection of the most basal electrode occurred at the beginning of the localization process. Furthermore, some electrodes were detected twice. In the present study, the localization accuracy of the method was investigated with different EAs. The decrease in the detection rate for the apical electrodes can be explained firstly by the nature of the pathfinding method used by Noble and Dawant [8]. The method adds a new node to the path at each iteration in the pathfinding process. This new node must satisfy certain constraints and is expected to minimize the cost function. However, some nodes are added to the path that do not match the desired position of a CI electrode. If a node is detected incorrectly, there is a high probability that subsequent nodes will also be detected incorrectly. Thus, the probability of correctly detecting an electrode at the beginning is higher than at the end of the search path, procedure last nodes in the path correspond to the apical electrodes. Furthermore, in our dataset it was found that where the often no candidate points are placed on the apical electrodes, which is a prerequisite to localize them (see Fig. 6, left). This effect was particularly clear for MED-EL EAs, because the apical electrodes have lower intensity in the CBCT image due to their design.

Generated candidate points for an example of the MED-EL’s electrode array Flex 28 for the method according to Noble and Dawant [8] with a threshold α₁ of 0.08% (right)

In our implementation of the algorithm by Noble and Dawant [8], the extraction of the centerlines from the binary mask was performed differently. It was analyzed whether it is possible to generate at least one candidate point on all electrodes with a threshold value proposed by Noble from the cumulative histogram. The threshold was set to α₁ = 0.08% as in Noble and Dawant [8]. As observed in Fig. 6 (right), no candidate points were generated on the apical electrodes. Moreover, it can be seen in Fig. 6 (right) that no artifact can be seen on CI electrodes 8,9,11 and 12. This means that the threshold derived from the cumulative histogram caused removal of relevant information from the ROI required for the localization of electrodes in EAs of MED-EL.

Therefore, it was necessary to re-parameterize the threshold from the cumulative histogram for each EA type. The aim was to select the threshold that maximized the likelihood that each CI electrode gets assigned at least one candidate point.

To ensure this, the manually determined electrode positions are used for the determination. To parameterize the threshold α₁, a range of values from 0.01 to 3% with a step size of 0.01% was investigated. This range was chosen so that there was a threshold for which all CI electrodes were included in the ROI and there were hardly any disturbing artifacts in the ROI. For each threshold, Fig. 7 shows how many electrodes have at least one candidate point within the radius of the GT position, where the radius is 0.9 mm (3 pixels). Table 4 shows the thresholds for each EA and the averaged thresholds across EAs.

Percentage of electrodes with at least one centerline candidate is allocated for each electrode array type. The right most subplot presents the average across each type of electrode array

Advanced cost function

For large insertion angles, the apical electrodes are close to the basal electrodes with a distance that is in the same range as the electrode spacing of the neighboring electrodes or inter-electrode distance. Paths that cross the EA from the apical to the basal side or vice versa might have a lower cost than pathways that arrange the electrodes in the correct order. In addition, Hachmann and Nogueira [7] reported that in the procedure of Noble and Dawant [8], the most basal electrode is often detected twice. Therefore, we proposed an extended cost function to prevent this behavior as follows:

C_R (redetection of CI electrodes) is the new extended cost term. To determine this cost term C_R, the Euclidean distances from a candidate point c to each node already found in the path p are first calculated. Then, the smallest difference diff_min between the candidate point c and a path node p_i is determined. This could be expressed mathematically as \({\text{diff}}_{{\text{min}}} = {\mathop {\min }\limits_{\,}} ({\text{diff}})\). Now diff_min is compared to the defined threshold value α₁₂, which describes the limit range for which no further child node c will be added to the path around a node found in the path. The threshold α₁₂ is set to 2/3 of the electrode spacing defined in the EA specification. If diff_min is less than α₁₂, a cost of infinity is assigned to C_R, and if diff_min is greater than this threshold a cost of 0 is assigned. The new cost factor C_R is integrated as a summand into the existing cost function of Noble and Dawant [8] as shown in Eq. (14).

Analysis of errors made by the algorithm to automatically localize cochlear implant electrodes

In the first step of the analysis, the results of the automatic localization algorithm were manually analyzed. In this way, errors that occurred during the automatic localization process can be classified. Table 5 shows these errors in percentage.

Three types of errors were considered: first, the electrodes can be detected twice. Second, the apical and basal electrodes can be confused during the detection process, in this case two subtypes of errors can be identified. In the first subtype, an apical–basal (a–b) confusion can occur (see Fig. 1, left). In this case, first basal electrodes are detected correctly, but afterwards it may happen that one of the basal electrodes is detected again. In the second subtype, a basal–apical confusion is considered (see Fig. 1, right). Here, an apical electrode is detected after the first basal electrode and the electrodes are then detected backwards. Third, nodes are incorrectly placed at locations corresponding with artifacts caused by, e.g., bones or wire leads in the localization procedure. Using the generalized threshold algorithm the third type of errors occurs in 55% of the CBCT data sets. This can be explained by the high threshold value used in the generalized threshold algorithm causing that many candidate points are located on artifacts and not on the actual CI electrodes. With the individualized threshold, this error occurs only in 12% of the CBCT data sets. The error is most pronounced with the Flex 16 EAs from MED-EL (67%) and with the Hybrid-L EAs from Cochlear (47%). For both types of EAs, the threshold value was chosen very high compared to the other EAs, generating many candidate points on artifacts. Thus, for this EAs it is more likely that nodes in the path are placed on artifacts and therefore detected as electrodes by the localization procedure. In addition, the Flex 20 EAs from MED-EL, Nucleus CI24RE and Nucleus CI624 also occasionally detected artifacts as nodes in the path. This indicates that the selected threshold also generates candidate points on artifacts. By extending the cost function, the error could not be minimized. This was expected, since the newly proposed advanced cost function does not aim at minimizing the detection of artifacts. However, for some data sets, the number of cases where nodes in the path were located on artifacts increased by 7 percentage points for MED-EL’s Flex 24, by 7 percentage points for MED-EL’s Flex 20, and by 27 percentage points for Cochlear’s Hybrid-L. When comparing the CBCT data sets in which artifact detection occurred and after extending the advanced cost function based on the individual threshold results, it was found that the advanced cost function eliminated localization errors related to double electrode detection, but resulted in artifact detection. On average, localization errors occurred for the generalized threshold and the individualized threshold algorithms for a similar number of data sets. Using the newly proposed advanced cost function, the detection error could be reduced. Errors occurred only in 3 CBCT data sets of the EA Hybrid-L from Cochlear. For these datasets, no improvement in localization could be achieved using the new advanced cost function. Compared to the generalized threshold algorithm, there was an increased confusion of apical–basal electrodes than with the individualized threshold. Using the advanced cost function, the apical–basal confusion were completely eliminated. When considering the data sets where basal apical electrode confusions were observed, no double detection of electrodes occurred. Thus, duplicate detection of electrodes was eliminated by extending the cost function in 98% of the CBCT data sets.

Statistical analysis

All statistical analyses were performed using the Friedman’s test and Dunn–Bonferroni post hoc tests. A p-value of less than 0.001 was chosen as the level of significance.