Automatic Localization of Deep Stimulation Electrodes Using Trajectory-based Segmentation Approach

TL;DR

This paper presents an automated, threshold-based segmentation method for accurately localizing deep brain stimulation electrodes in CT images, aiding precise VTA identification in Parkinson's disease patients.

Contribution

It introduces an adaptive threshold segmentation approach for robust electrode localization in noisy CT scans, improving accuracy over manual methods.

Findings

High noise tolerance demonstrated in CT image segmentation

Automatic threshold detection improves localization accuracy

Method facilitates better VTA estimation in DBS surgery

Abstract

Parkinson's disease (PD) is a degenerative condition of the nervous system, which manifests itself primarily as muscle stiffness, hypokinesia, bradykinesia, and tremor. In patients suffering from advanced stages of PD, Deep Brain Stimulation neurosurgery (DBS) is the best alternative to medical treatment, especially when they become tolerant to the drugs. This surgery produces a neuronal activity, a result from electrical stimulation, whose quantification is known as Volume of Tissue Activated (VTA). To locate correctly the VTA in the cerebral volume space, one should be aware exactly the location of the tip of the DBS electrodes, as well as their spatial projection. In this paper, we automatically locate DBS electrodes using a threshold-based medical imaging segmentation methodology, determining the optimal value of this threshold adaptively. The proposed methodology allows the localization of DBS electrodes in Computed Tomography (CT) images, with high noise tolerance, using automatic threshold detection methods.

Tables

Table 1. Table 1: Dimensions of the electrode trajectory in MNI space. These distances are calculated on the axial axis, taking into account the slice where the electrode tip is located and the slice where the DBS electrodes touch the patient’s skull.

Patient	Electrode	Electrode	Absolute	Volume size
ID	start (mm)	end (mm)	distance (mm)	axial (mm)
P1	-36	48.5	84.5	114
P2	-11	51.5	62.5	198
P3	-22.5	43	65.5	168
P4	-50	40	90	150
P5	-41.5	20.5	62	162
P6	2.5	71	73.5	156
P7	-12.5	66	78.5	221
Mean Vale	-24.429	48.643	73.786	167

Table 2. Table 2: Mean Squared Error (MSE) between Ground Truth and the coordinates of electrode trajectory generated by LEAD and proposed Method. Last, the average value and standard deviation for each method are reported for both electrodes.

Patient ID	Right Electrode		Left Electrode
	MSE OWN (mm)	MSE LEAD (mm)	MSE OWN (mm)	MSE LEAD
P1	0.049	81.68	0.35	54.38
P2	0.05	0.442	0.038	2.049
P3	0.04	0.91	0.1	0.83
P4	0.054	207.49	0.174	229.91
P6	0.24	3.13	0.081	3.93
P7	0.051	4.91	0.338	4.84
P8	0.448	29.307	0.197	125.12
Mean $\pm$ SD	0.134 $\pm$ 0.155	46.838 $\pm$ 76.66	0.182 $\pm$ 0.122	60.151 $\pm$ 87.737

Equations

$$Y=\left\{x_r\in R^{s\times a}:r\in 1,\cdots ,a\right\}.$$

$$\hat{\mathbf{Y}}=\mathbf{Y}\circ \mathbf{Z}.$$

$$\mathbf{B_{r}}=\left\{\begin{array}[]{ll}1&\mathrm{if\ }\mathbf{Y(x)}\circ\mathbf{Z(x)}\geq T_{r}\\ 0,&otherwise.\end{array}\right.$$

$$T=[K\cdot \sigma ]+\mu .$$

Full text

In this paper, we automatically locate DBS electrodes using a threshold-based medical imaging segmentation methodology, determining the optimal value of this threshold adaptively. The proposed methodology allows the localization of DBS electrodes in Computed Tomography (CT) images, with high noise tolerance, using automatic threshold detection methods.

1 Introduction

Parkinson’s disease (PD) is a degenerative condition of the nervous system, which manifests itself primarily as muscle stiffness, hypokinesia, bradykinesia, and tremor. In patients suffering from advanced stages of PD, Deep Brain Stimulation (DBS) neurosurgery is the best alternative to medical treatment, especially when they become tolerant to the drugs. DBS as a surgical treatment for PD symptoms was approved in 1998 by the United States Food and Drug Administration. This surgery involves the implantation of a medical device which consists of three components: Implanted Pulse Generator (IPG), DBS electrodes, containing one to four platinum iridium contacts, and electric conduction cables. These three elements are implanted surgically into the human body [1].

Research on DBS is currently being carried out in various centers around the world, and there is an extensive published scientific literature on the subject [2], [3]. DBS surgery produces a neuronal activity as a result of electrical stimulation generated by the electrodes. The quantification of neuronal activity is known as Volume of Tissue Activated (VTA) and can be interpreted as the amount of brain tissue that presents excitation or electrical response to stimulation generated by DBS electrodes. This VTA is a function of the DBS electrodes location within the brain. Therefore, the information regarding the position of these electrodes is used for the tuning of the stimulation.

Studies in patients with PD treated with DBS surgery show that the precision in the location of the electrodes within the desired anatomical region is directly correlated with motor enhancement during DBS [4], [5]. Although other factors influence the clinical results of DBS [6], the wrong location of the contacts is considered the most common cause of poor clinical response [7]. Consequently, knowledge about the location of the contacts becomes critical information in clinical routines to evaluate the effects of DBS and may help to define the best contact of the electrode that should be stimulated to treat chronic DBS, especially in the majority of current electrode designs, which have multiple stimulation contacts.

The location of electrodes used in DBS is still a challenging procedure, due to the presence of metal artifacts in the electrode and system cables for DBS [8], [9]. Dykstra et al. [10] worked on approach with DBS electrodes, performing contact coordinates extraction by visual inspection of post-surgical CT. However, this method is limited to a reduced number of brain areas, due to the time required for visual extraction. Lalys et al. [11] worked on identifying optimal sites for DBS, constructing an atlas from post-surgical images. For the segmentation of electrodes used in the DBS, they applied a fixed intensity threshold, whose value is not specified, to extract the lower voxel of the patient electrodes. Lalys obtained an error of 1.31 mm while finding the barycenter of the voxel extracted over a 10 mm extension, which then allowed them to get the location of the contacts using the geometric constraints of Medtronic’s 3389 electrode. However, they did not use actual electrode images for this location, but rather modeled the electrode axis within the atlas, using regression methods.

Motevakel and Medvedev [12] used a technique called “connect method”, which takes advantage of the linear shape of the electrodes and light beams reflecting the metal elements, resulting in a maximum error of 2.5 mm in the location of the electrodes used in DBS. However, their technique required that the two electrodes were present, which does not make them useful in the clinical cases where only one DBS electrode is introduced. To solve this limitation, some authors perform an artifact analysis for the electrode, based on semiautomatic algorithms [13]. Silva et al. [9] used a cloud-based tool [14] to estimate automatically the position of the DBS electrodes compared with anatomical structures, using an estimation of the rectilinear trajectory of the DBS electrode, based on the upper and lower tip. For the location of each cranial electrode, they used the distances in mm from each contact with the bottom tip, specified by the manufacturer. Horn et al. [7] determine the location of the electrode, through a semi-automatic process. They develop a Toolbox in MATLAB, which it has a precision range in the path location of two mm. Finally, they find the electrode using a straight-line model.

The medical image segmentation is a field of study that has been active for the last few decades, and results exist where the threshold varies proportionally to the average of the intensity of the CT [7]. The methodology here proposed to determine intrinsic parameters of each CT potentially allow the segmentation to be adjusted to each patient and to have resulted not so generic but adjusted to the noise, distortion and other variables that vary between each CT.

In this paper, we propose a methodology for automatic segmentation of the DBS electrodes, which adjust to the intensity and rotation parameters of each CT scan. This method includes the CT image register with a clearly delimited anatomical atlas, which allows the creation of a Region of Interest (ROI) fully adjustable to the anatomy of each patient, followed by an intensity value that acts as a threshold for the segmentation, which is performed in this ROI. The performance of this methodology is measured by the comparison with a Ground Truth, using a Euclidean distance.

2 Materials and Methods

This section presents the materials and methods used for the development of the proposed methodology. The Figure 1 shows the general diagram, where the gray boxes represent the medical images as inputs, blue boxes represent state of the art techniques proposed by others authors, and green boxes represent the new contributions made it. Finally, the red box represents the centroids defining the trajectory of DBS Electrodes, as main output from the method.

2.1 Atlas Originated

Due to significant variations in orientation in different CT images, which will depend on from the CT scanner in which they are performed, a single ROI cannot be generated for all patients. Therefore, this research proposes an iterative and autonomous method, where an ATLAS is created for a patient, defining a region of optimal and precise interest for that patient, determined by well-defined anatomical coordinates. This ROI should cover most of the DBS electrodes trajectory but should minimize the presence of bone structures within it.

initially, to determine a geometric approximation of this ROI, the dimensions of the electrode trajectory in MNI space are taken into account, as shown in Table 1. This lists the dimensions of the electrode trajectory in the patients analyzed, so it provides significant geometrical constraints to generate the ROI defined by the coregister with Atlas.

2.2 Coregister Atlas and CT MNI

For this coregister, which represents the third system of Figure 1, the software Elastix [15] is used. The MRI image registration process is performed by a linear process with a rigid register, and then a register based on an affine transformation.

Let $\Psi=\left\{x_{r}\in\mathbb{R}:\mathbf{r}=(i,j,k)\right\}$ be an structured set of intensity voxels of image CT normalized to MNI space, where $x_{r}$ is the Hounsfield intensity measured at location $\mathbf{r}\in\mathbf{R}\subset\mathbb{N}^{3}$.

[TABLE]

Where $r$ is the number of slices will be analyzed. Let $\mathbf{\hat{Y}}$ be an ROI, obtained by multiply;

[TABLE]

Where the function $\mathbf{Z}$ is an ROI defined as $\mathbf{Z}\in[0,1]^{M\times N\times O}$.Finally, we obtained a binary image $\mathbf{B}$, defined as:

[TABLE]

[TABLE]

Where $T_{r}$, is the threshold calculated for each slice in the axial axis.

2.3 Adaptive threshold

The electrode trajectory was identified using threshold segmentation. For this, we used the approximation by Rodriguez et al. [16] and Otsu et al. [17], where they describe the Local Adaptive Thresholding method. Because there exist a significant level of variation in the intensity of the gray-scale level within the range of $0-4000$ HU, the global threshold does not work well here. Local threshold works best for this case where the image is divided into slices, so we find the single threshold for each slice. The threshold value will depend on regional statistics such as variance and mean of the image. We get the threshold $T$ by adding to the intensity average $\mu$ of the slice in the CT image, the amount $K\times\sigma$ of the Eq. 5:

[TABLE]

Where $\sigma$ and $\mu$ are, respectively, the standard deviation and mean of intensities values of pixels, within the local window. $K$ is a variable that can be determined to obtain the best segmentation. However, approaching the threshold in this way, with the proposed value for $K=0.9$ by [7], it can be shown that all or most of the voxel intensity values of DBS electrodes can be calculated below the threshold, which will not allow a correct segmentation. To find the optimal value of K parameter that provides the best segmentation, we next performed an exploratory analysis, where K value is varied by of $\textpm 0.01$, in the range $[0.6-9.5]$. After was evaluating the Euclidean distance that defines the segmentation for each of these values of K. The upper limit for the range is determined in $9.5$; since from 9.5, the error initiates an upward trend, and this is a fact that is repeated in patients of the database.

2.4 Validation

Because it is not desired to compare the volume of the electrode, but only its centroid, we use Euclidean distance, which establishes the separation distance between the centroid of Ground Truth, and the centroid calculated with the proposed methodology. To compare the results with the results obtained in [7], we must decrease the resolution of the planned trajectory. So, the trajectory with a resolution of $0.5$ mm is again sampled so that in the end it has a resolution of $0.63$ mm (LEAD). For this step, we used the concept of linear interpolation.

2.5 Dataset

The dataset was collected at the Instituto de Epilepsia y Parkinson del Eje and Instituto Neurológico de Colombia. It is composed of eight (8) computed tomography scans, applied to an equal number of patients, all diagnosed by a specialist as patients with Parkinson’s Disease. These imaging studies by Tomography have an average resolution of 1.45 mm. In order to get the Ground Truth, the centroids of the electrodes are manually labeled in the slices containing electrodes. In this study, were considered more than 1600 by hand labeled centroids, belonging to slices that contain electrodes.

3 Results and Discussion

In Figure 2 is shown the mean Euclidean distance (ED) that defines the segmentation accuracy for six patients, using the proposed method for K values from the range [$0.6-9.8$]. The ED presented by LEAD, on average, for the six patients, is shown with red stars. This distance for LEAD is calculated only at K = 0.9. However, an extension to this line is made, using red stars, to show that the ED that generates the method, is consistently inferior to the ED of LEAD.

In more than 50 $\%$ patients, the methodology has a lower Euclidean distance, even with the same K value as the one used in LEAD, then showing the improvement in the precision provided by the proposed method. The unique conditions as noise and intensity for each patient made that a single K value cannot be defined as the only appropriate value, which will be able to use for all patients, we found that a qualitative analysis of the range should be performed. In this qualitative study, a range of optimal K values is determined, for which the Euclidean distance is minimized for the database patients, thus finding the best possible general segmentation.

From Figure 2 we can deduce that all K values higher than four, have a lower Euclidean distance than LEAD. However, to make more precise the location of the K range, it is placed as a condition that all K values within this range should allow a segmentation whose Euclidean distance with Ground Truth is less than $1$ mm for all patients. Finally, we obtain that the range of values for K selected for the methodology is [$7.27-8.5$].

3.1 Proposed Atlas for coregister

Due to significant variations in orientation in different CT images, which will depend on from the CT scanner in which they are performed, a single ROI cannot be generated for all patients. Therefore, this research proposes an iterative and autonomous method, where an ATLAS is created for a patient, generating a region of optimal and precise interest for that patient, determined by well-defined anatomical coordinates. This ROI should cover most of the DBS electrodes trajectory but should minimize the presence of bone structures within it.

Then the Atlas for the ROI coregister is defined anatomically, precisely and easily reproducible, as follows:

In the sagittal plane, there are no variations; the only limit is the range [$-40,+40$] mm. 2. 2.

In the caudal part, the axial axis has values of [-49, -39.5] mm, while the coronal axis is limited between [-30, +3] mm. 3. 3.

In the immediately upper part, the axial axis is limited in the range [-39, +1.5] mm and the coronal axis is limited between [-30, +22] mm. 4. 4.

In the most cranial part of the Atlas, the axial axis is limited in the range [$+2,+31$] mm and the coronal axis is limited between [-30, +32] mm.

This ROI is 80 mm deep in axial axis, has a width of 80 mm in sagittal axis, and in coronal axis, the length changes because the electrode trajectory travels mainly in this axis. It was set up from Patient 4 CT, this because it is the patient who presents the greatest depth of the electrode in the axial axis. So, it is the patient who could introduce more significant measures of error and presence of cranial structures. In Figure 3 the resulting mha format image (green) is shown after the coregister between ATLAS and MNI P3 CT. It is observed that the atlas is well delimited, minimizing the cranial area in analysis, and allowing to cover most of the electrode trajectory.

In table 2 is shown the Mean Squared Error (MSE) between the trajectory coordinates of both electrodes, delivered by the proposed method and the LEAD method, and the Ground Truth coordinates. Last, the average and standard deviation of this error are reported. It is observed that the proposed method consistently has a reduced error compared to the error generated by the LEAD method. It should be noted this may be due to the LEAD method is not an automatic method, so a manual correction must be made in the height on the axial axis by the specialist. In the error measure presented here, there is no kind of height correction was performed, which explains the presence of an MSE of up to 207 mm in some patients for the LEAD method.

In Figure 4 is shown the 3D reconstruction of DBS electrode trajectory provided by the proposed and LEAD method placed over the image of the electrodes (blank). It is observed that the LEAD method (red) does not define well the electrode trajectory. The Figure 4(a) shows the best result obtained by the LEAD Method in the electrode segmentation and the figure 4(b) the worst case is shown. It is observed that although LEAD follows the trajectory correctly in the coronal and sagittal axes, it does not have concordance with the trajectory along the axial axis. This lack of accuracy of LEAD is because it is not an automatic method, so it requires a specialist adjustment in the axial axis depth. Instead, the proposed method segments the electrode trajectory automatically (green), without supervision or correction of any kind. The electrode trajectory obtained by applying this methodology matches with an accuracy of $0.2$ mm the right trajectory of the electrodes, thus preliminary validating our research.

The Figure 5 shows the time comparison of the segmentation process of the DBS electrodes, both in Toolbox LEAD as using the methodology proposed here. It is demonstrated that the proposed method takes around $28\%$ of the time it takes LEAD to arrive at a result. That is, based on the computation times in Fig. 5, there is an improvement of $78$% in the computational time taken by the method on the time it takes the Toolbox LEAD. The computational times were measured in a dual-core i5 2450-M 2.5 GHz processor, based on the Sandy Bridge architecture, with 12 GB of RAM, Intel HD Graphics 3000 integrated video card and 128 GB internal SSD.

The accuracy DBS electrodes segmentation is an important task, particularly since a spatial and geometric reference point for any tool that intends to generate information of physical or chemical processes involving these DBS electrodes. That is why this paper provides a new methodology that allows DBS electrodes segmentation to be performed with an accuracy of less than $0.2$ mm, in the eight patients in the database. Due to a limited number of patients for to validate the results, it is important to corroborate these results with larger databases, to verify the proper functioning of electrodes segmentation. Regarding state of the art, the present work makes a contribution by proposing a methodology that has proven to perform the segmentation of the DBS electrodes automatically and accurately. This advanced method suppresses the independence of specialist that have other tools since the electrodes trajectory are not required manual corrections in any of three anatomical axes.

4 Conclusions

In conclusion, we provide here a methodology of segmentation of medical images for automatic localization of electrodes used for the surgery of Deep Brain Stimulation, a procedure that is performed in the treatment of Parkinson’s disease. The methodology proposed has been shown to adapt to the patient’s conditions regarding noise, distortion and other parameters of the CT. Results indicate that the method can locate with an accuracy of less than 1.5 mm per slice, the effective path followed by the DBS electrodes, improving what is found in state of the art. The proposed method has positive clinical implications, such as that specialist should not provide feedback in the process of the DBS electrodes segmentation, since the method proposed here in the dataset used has proved to be automatic and accurate. Additionally, the computation time of DBS electrodes segmentation has dramatically decreased, reaching reductions of up to $78\%$ in calculation time of the trajectory, compared with previous approaches.

Future work can be focused on:

• •

Future work 1. Validate the proposed methodology with a larger number of patients, with the objective of patenting both the methodology and the software, for clinical use.

• •

Future work 2. To implement the necessary changes in the source code that allow a parallel execution, using concurrent programming and by threads; all this to improve the execution time.

• •

Future work 3. To implement the methodology proposed in free software, such as the ITK library, following a pipelining method, which allows executing all the computer libraries from a single software environment.

• •

Future work 4:

Acknowledgments. This work was developed within the project 1110-657-40687 “Estimación de los parámetros de neuro modulación con terapia de estimulación cerebral profunda en pacientes con enfermedad de Parkinson a partir del volumen de tejido activo planeado”, with financial support of Colciencias (Colombia). Authors were also supported by the 706 agreement, “Jóvenes Investigadores e Innovadores”, with the proposal “Desarrollo de una metodología para localización automática de electrodos utilizados en DBS (estimulación cerebral profunda) para el tratamiento de la enfermedad de Parkinson”, funded by Colciencias (Colombia). Finally, the authors are thankful to the ‘Grupo de Investigación en Automática’ ascribed to the engineering program at the Universidad Tecnológica de Pereira, and M.Sc. H.F. García for technical support.