Automatic segmentation of stereoelectroencephalography (SEEG) electrodes post-implantation considering bending

Automatic segmentation of stereoelectroencephalography (SEEG) electrodes post-implantation... Purpose The accurate and automatic localisation of SEEG electrodes is crucial for determining the location of epileptic seizure onset. We propose an algorithm for the automatic segmentation of electrode bolts and contacts that accounts for electrode bending in relation to regional brain anatomy. Methods Co-registered post-implantation CT, pre-implantation MRI, and brain parcellation images are used to create regions of interest to automatically segment bolts and contacts. Contact search strategy is based on the direction of the bolt with distance and angle constraints, in addition to post-processing steps that assign remaining contacts and predict contact position. We measured the accuracy of contact position, bolt angle, and anatomical region at the tip of the electrode in 23 post-SEEG cases comprising two different surgical approaches when placing a guiding stylet close to and far from target point. Local and global bending are computed when modelling electrodes as elastic rods. Results Our approach executed on average in 36.17 s with a sensitivity of 98.81% and a positive predictive value (PPV) of 95.01%. Compared to manual segmentation, the position of contacts had a mean absolute error of 0.38 mm and the mean bolt angle difference of 0.59 resulted in a mean displacement error of 0.68 mm at the tip of the electrode. Anatomical regions at the tip of the electrode were in strong concordance with those selected manually by neurosurgeons, ICC (3, k) = 0.76, with average distance between regions of 0.82 mm when in disagreement. Our approach performed equally in two surgical approaches regardless of the amount of electrode bending. Conclusion We present a method robust to electrode bending that can accurately segment contact positions and bolt orien- tation. The techniques presented in this paper will allow further characterisation of bending within different brain regions. Keywords Epilepsy · SEEG · Automatic segmentation · Bending Introduction Electronic supplementary material The online version of this article (https://doi.org/10.1007/s11548-018-1740-8) contains supplementary Epilepsy is a disease characterised by an enduring predis- material, which is available to authorized users. position to generate epileptic seizures and affects 1% of the population [8]. A third of patients develop chronic refractory B Alejandro Granados focal epilepsy and neurosurgery may provide a cure [9]. alejandro.granados@ucl.ac.uk Brain imaging is fundamental in a typical neurosurgi- Wellcome/EPSRC Centre for Interventional and Surgical cal evaluation for determining the epileptogenic zone (EZ) Sciences, UCL, London, UK National Hospital for Neurology and Neurosurgery, London, The First Affiliated Hospital of Xian Jiaotong University, UK Xian, People’s Republic of China Department of Clinical and Experimental Epilepsy, Institute Vickie and Jack Farber Inst for Neuroscience, Thomas of Neurology, National Hospital for Neurology and Jefferson University, Philadelphia, USA Neurosurgery, London, UK Dementia Research Centre, Department of Neurodegenerative Department of Statistical Science, University College Disease, UCL Institute of Neurology, London, UK London, London, UK 123 936 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 with modalities including structural and functional MRI trodes are modelled fitting a polynomial function and then (e.g. T1/T2-w, FLAIR) and PET [9]. If the EZ is not aligned to a common coordinate system reporting mean identifiable, invasive electroencephalography (EEG) record- deviation that varies from 0.92 to 2.0 mm. However, they ings are performed in the form of stereo-EEG (SEEG) have mostly focused on fitting trajectories using polynomi- or subdural grid insertion. SEEG is a procedure in which als rather than computing the amount of electrode bending multiple electrodes are stereotactically inserted to identify and have not considered the reasons of bending within the the seizure onset zone [21]. Accurate placement of elec- brain anatomy. Although Lalys et al. [13] looked at the rea- trode contacts is important for safety, interpretation of the sons of bending (mainly due to brain shift) by computing recorded electrical signals, and subsequent resection plan- a local and mean curvature index over the entire length of ning [21]. Planning of electrode implantation is crucial for DBS electrodes, the index provides no information about the avoiding blood vessel damage and subsequent intracranial direction of bending. Unlike SEEG procedures, where 8–14 haemorrhage (which occurs in 1–2% of patients), and auto- electrodes are inserted, DBS electrodes are typically inserted matic computer-assisted multiple trajectory planning tools bilaterally and the contacts are very close to the tip. To dis- have been proposed [17,18]. However, intraoperatively, entry criminate between contacts located in white or grey matter, point (EP) accuracy can be affected by misregistration of the Arnulfo et al. [1,16] compute the distance from each contact neuronavigation system, inaccurate alignment, and deflec- to grey–white matter interface. tion during drilling, whereas target point (TP) errors may be caused by the angle at which the electrode passes through skull, deflection of the electrode at the dura or within the Contribution of this paper brain, the rigidity of the electrode, and the depth to which a guiding stylet is inserted [3,21]. Robotic systems have Our main motivation is to automatically segment SEEG con- been introduced to improve EP implantation accuracy [3,7]. tacts and bolts (Ad-Tech Med Instr Corp, USA) relative to However, TP displacement is the main source of error and the anatomy whilst accounting for electrode bending along understanding why and how electrodes bend may help pre- its trajectory at contact positions rather than as a result of TP dict final TP positions during surgical planning and improve displacement. Our algorithm (Fig. 1) allows estimating not EZ localisation [22]. only the position of contacts but also the direction of the bolts Furthermore, it is convenient to have a rapid and reliable inserted into the skull since the angle of the bolt with respect scheme for segmenting contacts, assigning their anatomical to the scalp surface normal is a measure of post-implantation location when interpreting SEEG studies and for guid- accuracy. We quantify local and global bending by means of ing definitive surgical resections. Automatic segmentation electrodes modelled as elastic rods in position-based dynam- approaches have been proposed for SEEG [2,14,16] and deep ics and validate our methods in 23 post-SEEG cases (224 brain stimulation (DBS) [5,10,11] implantation. Arnulfo electrodes, 1843 contacts) comprising two different surgical et al. [2] used post- implantation CT (threshold = 1600) approaches (placing a guiding stylet close to or far from the co-registered with MRI to segment electrodes based on a TP). geometrical-constrained search. They randomly generated different scenarios for 1–15-mm displaced TP in an experi- mental study and reported accuracy of 10% of false negatives Methods (FN) and 7% of false positives (FP) for a maximum dis- placement of 15 mm. However, bending may occur at any Input images point along the electrode’s trajectory. Meesters et al. [14] co-registered the CT (threshold = 500 HU) with MRI and A post-SEEG implantation resampled CT and an MRI T1 extracted guiding screws with a multi-scale filter whilst deter- images are rigidly co-registered using NiftyReg (v1.5.43) mining likely tip locations within a wedge-shape region. [15]. From the MRI image, we obtain the parcellation of However, manual adjustments took between 10 s and several brain anatomy via NiftyWeb (GIF v3.0) (Fig. 2)[4]. minutes, and reported deviations of the tip and their method did not account for electrodes bending. Additionally, these Identification of anatomical masks methods relied on pre-operative plans and were tested only on one electrode type. We use the MRI and the parcellation to create regions of Hubsch et al. [10,11] proposed an automated algorithm interest that are used to identify contacts, bolt heads, and the reconstructing full electrode trajectory whilst accounting for section of the bolt crossing the scalp/skull, which we refer DBS electrode bending from CT scans. A convex hull brain mask is extracted using thresholds, and the largest connected https://github.com/InteractiveComputerGraphics/PositionBasedDyn components are skeletonised [10]. Trajectories of 11 elec- amics. 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 937 Fig. 1 Flow chart of algorithm pipeline Fig. 2 Input images: a post-SEEG implantation CT, b MRI T1, and c parcellation Fig. 3 a Axial, b sagittal and c coronal planes showing computed masks of the brain (cyan), skull and scalp (yellow) together with the result of connected components filters of contacts (red), bolt head (green), and the section of the bolt crossing the skull (blue) as bolt body.First,a BinaryThresholdImageFilter A BinaryThresholdImageFilter is applied to the is applied to the parcellation to create a mask of intracra- MRI to create a mask of the scalp BI with a lower scalp nial space BI , i.e. with a threshold t in the range of threshold equal to t . We use morphological operators brain brain scalp 4 ≤ t ≤ 208. We apply a method similar to Dogdas et al. to combine BI and BI and apply a closing fil- brain brain scalp [6], which we describe herein for completeness. We compute ter with a ball structuring element (radius = 10) to obtain a a skull threshold t from the MRI as the mean of the inten- mask of the head, i.e. BI = (BI ∪ BI )  B , skull head scalp brain 10 sities of the nonzero voxels that are not brain as an empirical and a mask of the skull, i.e. BI = BI ⊕ BI , skull head brain measure to split the low- and high-intensity regions, followed after applying an XOR morphological operator on the result by a scalp threshold t as the mean of the non-brain vox- (Fig. 3). scalp els above the skull threshold (∀I (x , y, z) ≥ t )to MRI MRI skull identify the transition between the head and the background. 123 938 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 Table 1 Geometrical analysis, Geometrical analysis Discriminant analysis μ(σ ), and discriminant analysis Number of Pixels Elongation Roundness Number of pixels Roundness of bolt heads and contacts Bolt head 329.4 (183.5) 2.51 (0.59) 0.63 (0.06) > 100 [0.4, 1.0] Contact 9.7 (6.6) 2.52 (1.27) 1.10 (0.06) [3, 50] Fig. 4 Search strategy given the direction of the bolt and constraints (distance and angle) Segmentation of electrode bolts and contacts (30 ) (Fig. 4), constraints which favour assigning contacts in the direction of the bolt during a first pass. Amask BI is created from a BinaryThreshold post ImageFilter applied to the post-op CT with lower thresh- Automatic segmentation of electrodes old t = (0.52) ∗ max (I (x , y, z)). BI is used to CT CT post identify full bolts (BI ) with at least a minimum of 200 bolt The main steps of our algorithm include: pixels. Three subsections are identified: the head of the bolt which is outside the patient’s head (BI ∩¬BI ), the 1. Initialisation All segmented contacts are initially labelled bolt head as ’available’ and stored in a pool. Given a bolt head body (BI ∩ BI ), i.e. section crossing the skull, and the bolt skull position (x ), the closest bolt bodies (8 ≤ x − x ≤ tip (BI ∩ BI ). Lastly, contacts are identified within h h b bolt brain 25 mm) and the closest contact ( x − x ≤ 50 mm) the brain whilst excluding bolt tips ((BI ∩ BI ) ⊕ h c post brain are identified in order to narrow the search down to only BI ). We applied a ConnectedComponentImage boltTip Filter to the masks and a LabelImageToShape those relevant. 2. Contact search strategy For each bolt head, the contact LabelMapFilter to the blobs to get their centroids and geometrical properties before conducting geometrical anal- search strategy is executed initially with the closest bolt body (1st pass search) and subsequently with alternative ysis to identify discriminants of segmentation (Table 1). We detected contacts with blobs that were within a range of num- bolt bodies if no contacts have been assigned. Although rare, bolt bodies may not be segmented and a direction of ber of pixels ([3, 50]) and bolt heads with blobs that had a minimum number of pixels (≥ 100) and were within a range search cannot be computed. Therefore, the contact search strategy is called again with the closest contact position of roundness values ([0.4, 1.0]). rather than a bolt body position. 3. Project remaining contacts in pool For electrodes con- Contact search strategy taining at least one contact, we compute the minimum distance between an available contact in the pool and Given a bolt head (x ) and its closest bolt body (x ) posi- a line formed by the positions of the bolt head and the h b x −x b h tions, we compute the direction of search (nˆ = ) electrode tip. The contact is assigned to the electrode if x −x b h and iteratively compute a number of points x given a maxi- and only if its distance to the closest point x to the line p p mum electrode length (90 mm) and a step size (1 mm) in the (tangent to the line) is below a constraint (5 mm) and x direction nˆ . An available contact x is assigned to the elec- remains along the line or in a position of the line 20% 0 c trode if and only if it is located below a distance constraint extended from the tip, i.e. within an interpolation range from x (5 mm) and the angle between the previous direction of [0.0, 1.2] to project contacts that are further from the nˆ and the current direction nˆ is below an angle constraint currently identified tip of the electrode. 0 c 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 939 Fig. 5 Modelling of electrodes as elastic rods. Bolt head (green) and body (blue) with contacts (red) modelled as point particles and ghost particles (cyan) created orthogonally along the electrode with material frames located between contacts 4. Predict contacts in the bolt region For a given electrode, electrodes (N = 109 contacts), (b) manually identify the we compute the most common segment along the elec- tip and head of the bolt of a random subset of electrodes trode based on the distances between subsequent contacts (N = 95 bolts), (c) confirm the correct number and location rounded to the closest integer. Based on electrode spec- of contacts and electrodes (N = 23 cases), and (d) identify ification, we infer the type of electrode depending on the TP anatomical region (N = 222 electrodes). the order of the segments and specify contact spacing. We then compute the direction from the last contact x towards the bolt head x and create new contacts up to 21 Results mm before the bolt head position to segment only those contacts closer to the skull. Interface Bending estimation We implemented our algorithms in C++ using MITK and ITK as well as a GUI in Qt to allow clinicians to adjust the To quantify electrode bending, electrodes are modelled as automatic segmentation if needed (Fig. 6). On average, our elastic rods using the Cosserat model proposed by [19] method executed in 36.17 s (N = 23, σ = 15.7), faster than and then implemented by [12] in position-based dynam- manual segmentation. ics. Electrode contact positions are represented as linked particles with ghost particles located orthogonally half-way between contact pairs (Fig. 5). A material frame is created Performance between contacts with a unit vector (d = X − X ) 3 c c n−1 n aligned tangentially to its centreline followed by two addi- Of a total of 224 electrodes (1843 contacts), 29 contacts were tional orthonormal vectors, (d = d × (X − X )) and 2 3 c c n−1 segmented but not assigned to any electrode due to: (a) three ˆ ˆ (d = d × d ) chosen to lie in the principal direction of the 1 2 3 bolt heads that were not automatically segmented (17), (b) cross section. We compute the rate of change of two consec- no segmented contacts close to them (5), and (c) due to one utive frames, namely a Darboux vector , to describe local incorrectly assigned contact to a bolt head (7). On average, bending at the contact points [12,20]. Along the electrode, Tr P Tr P the sensitivity ( ∗ 100) and PPV ( ∗ 100) of Tr P+FN Tr P+FP values are then accumulated to quantify global bending. We our approach was μ = 98.81%; σ = 2.04 (false-negative then use the parcellation to report the region at which each rate of μ = 0.124; σ = 0.02) and μ = 95.01%; σ = 6.73 contact is located and report all those regions that the elec- (false- positive rate of μ = 0.059; σ = 0.09), respectively trode passes through. Lastly, contact displacement and depth (Fig. 7 bottom), finding no statistical significant difference are estimated with respect to a rigid electrode with position between data sets of the two surgical approaches, i.e. placing of contacts projected along the direction from the bolt head to a stylet far from or close to target point. To illustrate our the last contact (X ) at distances subject to electrode speci- results, Fig. 8 shows two cases correctly identified (a, b) fication. along two worst cases (c, d). Validation We asked two neurosurgeons and one clinical scientist to http://mitk.org. (a) manually segment the contacts of a random subset of http://itk.org/. 123 940 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 Fig. 6 Automatic segmentation of electrodes interface and GUI for manual adjustment Validation Computed contact positions, bolt angles, and regions of anatomy are compared with those manually segmented in a subset of cases (Table 2). – Contact position Compared to the manual segmentation done by a clinical scientist (M1) and a neurosurgeon (M2), we found that the contact location of our auto- matic segmentation approach had a mean absolute error (MAE) of 0.38 and 0.40 mm, respectively, and a root- mean-square deviation (RMSD) of 0.45. The distance of contact positions between both manual segmentations was on average μ = 0.37 mm (σ = 0.22). We found no statistical difference when comparing the distances from automatically computed contact position to those positions obtained via manual segmentation (paired dif- ferences: μ = 0.036, σ = 0.21). – Bolt angle We found that the angle of bolts between automatic and manual segmentation (Fig. 9) by M1 and ◦ ◦ M2 differed on average by 0.59 and 0.22 , respec- tively, with pair samples strongly and positively corre- lated (Pearson correlation) and with strong reliability (Cronbach’s alpha). We study the displacement error d = sin(θ ) ∗ l at the tip of the electrode error error e caused by this angle difference θ and the length error of the electrode l within the brain and define a maxi- mum tolerance value T = 2.29 mm related to contact length. On average, d at the tip of a rigid elec- error trode caused by the angle difference is μ = 0.68 mm Fig. 7 Top: Number of contacts initially segmented and assigned to for M1 and μ = 0.72 mm for M2. We found 3 out- electrodes via a bolt head and bolt bodies association (step 2), pro- jected (step 3), predicted (step 4), and left unassigned and available liers above T for M1 ([2.46, 4.79] mm) and 8 outliers in pool. Bottom: Number of contacts correctly segmented (TrP—true for M2 ([2.37, 5.48] mm). Given T , a non-inferiority positives), wrongly segmented (FP) and missed (FN) in 23 data sets test indicates that 0.68 mm is an estimate of d with error (order by number of electrodes) 95% of CI (0.431–0.926) after accounting for clustering 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 941 Fig. 8 Examples: a segmented bolt head and contacts of electrodes overlaying CT; b contacts predicted at the skull and scalp level; c 22 FP (red marks along the skull); and d our worst case with 3 contacts not segmented due to crossings and 4 FNs using a patient-level random effect. Figure 10 shows of electrodes (μ = 0.49; σ = 0.34) compared to the bend- an example with electrodes automatically segmented ing observed in electrodes that had a stylet placed close to TP and their corresponding rigid electrodes (lighter colours) (μ = 0.31; σ = 0.18), a difference which was statistically computed using the direction of bolts automatically significant, t (222) = 5.36, p < 0.01. segmented. – Regions of anatomy We also ran a intra-class correlation Generalisability and robustness two-way mixed effects model with average measures and found a strong agreement when identifying the anatomi- Three SEEG post-resection cases using SEEG DEPTH elec- cal region of the brain at the tip of the electrode between trodes (PMT Corp., USA) were obtained from the Vickie our algorithm and that done manually by two neuro- and Jack Farber Institute for Neuroscience (Thomas Jeffer- surgeons, ICC (3, k) = 0.76, p < 0.001. When in son University) to assess generalisability and robustness of disagreement, the average distance between regions was our algorithm. We observed a reduced average performance 0.82 mm (σ = 0.78), a distance below contact size. (sensitivity = 69.7% and PPV = 82.6%) due to the following Furthermore, it is estimated that electrode contacts elec- factors (Fig. 11a–c): (a) smaller bolt heads (our parameter of trically sample regions of grey matter within a 3 mm minimum number of pixels of bolt heads could be adjusted), radius. Any discrepancy in identified anatomical regions (b) contacts being very close to each other and merged as sin- below this is therefore not clinically significant. gle blobs (addressed by adopting optimal oblique resampling used for DBS electrodes [11]), and (c) electrodes inserted Bending estimation deeply (our parameter of maximum electrode length could be adjusted to account for this). Despite this, our algorithm In order to study whether Darboux vectors are a represen- was agnostic of electrode types and implantation planning tative measure of bending, we look into the relationship and was robust in post-resection data sets. (Pearson correlation) between global bending and the fol- We randomly chose 3 of our data sets to test the method lowing variables: accumulated displacement of contacts (r = proposed in [2,16] and implemented in SEEGA (Slicer 0.532, p < 0.001), length of electrode inside the brain tissue v4.6.2). We configured electrode types based on electrode (r = 0.373, p < 0.001), amount of white matter traversed specification, used the implantation plan (EP and TP) as fidu- by the electrode (r = 0.257, p < 0.001), and bolt angle cials and imported the CT image. We modified SEEGA to (r = 0.189, p = 0.045). Of the two surgical approaches, use the same threshold that we computed in our algorithm placing a stylet far from TP resulted in larger global bending for consistency and because the default threshold computed 123 942 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 Table 2 Validation between manual and automatic segmentation Measure Manual vs automatic Statistical test Result Contact position MAE (μ, σ , IQR) M1: 0.38 mm, 0.24, 0.22 Paired t test t (106) = − 1.756, p = 0.82 N = 109/1843 M2: 0.40 mm; 0.22, 0.26 Pearson correlation r = 0.454, p < 0.001 Cronbach’s alpha 0.615 MAE (x , y, z M1: (0.14, 0.15, 0.27) mm components) M2: (0.17, 0.15, 0.26) mm RMSD M1: 0.45 M2: 0.45 Bolt angle mean angle difference M1: 0.59 (1.27) Paired t test t (94) = − 4.54, p < 0.001 M1: N=95/224 M2: 0.22 (1.53) t (112) = 1.533, p = 0.128 Pearson correlation r = 0.991, p < 0.001 r = 0.985, p < 0.001 Cronbach’s alpha 0.995 0.992 M2: N = 113/224 displacement error at M1: 0.68 mm, 0.81, 0.83 Non-inferiority test CI = (0.431, 0.926) first contact due to tolerance = 2.29 mm M2: 0.72 mm, 0.84, 0.81 angle difference (μ, σ,IQR) Regions of anatomy region of anatomy at Intra-class correlation 0.76, p < 0.001 N = 222/222 first contact distance between 0.82 (0.78) mm regions when in disagreement Fig. 9 Bolt angles. a Bolt from post-CT image and b manual iden- fication (rigid electrode shown in blue) of bolt direction of an outlier tification of the direction along bolts by a clinical scientist. c, d case with angle difference of θ = 5.73 and displacement error at error (Inconspicuous) comparison of manual (pink) and automatic identi- the tip of d = 4.79 mm error by SEEGA resulted in segmentations errors. We observed noticed that most of the contact positions were not accurate an average sensitivity of 82.9% and PPV of 65.3% (97.3% (Fig. 11d), we only report performance on the number of con- and 98.2%, respectively, using our algorithm). Whilst we tacts that were incorrectly or not segmented. Similarly to our 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 943 Fig. 10 Anatomical regions traversal (top) and contact displacement (bottom) of automatically segmented electrodes with respect to a rigid electrode computed based on bolt direction approach tested with data from other centres, SEEGA might We found that FP were located in the inner surface of the perform better after fine-tuning parameters. Both algorithms skull and were caused by pixel size inaccuracies of CSF have parameters the user must select to guarantee optimal regions overlapping with bone structure. The performance performance. of our algorithm is similar to previous approaches although [2] only considers displacements at the tip of the electrode with no details of displacement of other contacts along the Discussion electrode and [14] uses a very small sample size assuming rigid electrodes. Automatic segmentation Compared to the search strategy by Arnulfo et al. [2], our algorithm uses a higher angle constraint (30 rather than The automatic segmentation of bolts is typically overlooked 10 ) because we use instead the bolt direction to search for in the literature and could be used to report accuracy errors contacts rather than a direction from previously segmented caused by differences in angle with respect to planning. contacts. It is also clinically relevant to accurately segment We use bold direction to search for contacts with neither the position of the contacts closer to the skull to ensure grey prior information of electrode type nor implantation plan- matter at the cortical entry is adequately sampled, but these ning. Compared to previous work [2,14], we use a factor might be difficult to segment in the bolt region. Therefore, of maximum intensity from CT images rather than a con- contacts are predicted in this region after inferring electrode stant. Regions of interests were used to segment position type. Further conditions would need to be included in this step of contacts and bolts based on geometrical properties with to support more electrodes from different manufacturers. To their centroids equivalent to the signal peaks found in [11] cope with electrodes crossing, we initially used geometri- and more generally in the literature. The choice of intensity cal features to identify large blobs that relate to more than threshold and constraints favour few incorrectly segmented one contact for splitting. However, the resulting position of contacts (FP) over missing contacts (FN), since these can contacts was not good enough to make the method fully auto- be easily discarded by surgeons during manual adjustment. matic, so we rely on manual adjustments that can be quickly 123 944 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 Fig. 11 Generalisability and robustness tests. Our proposed algorithm electrode with insertion depth of 110 mm); Our data in SEEGA: d seg- using data from a different centre: a with smaller bolt heads, b con- mented contact positions (green fiducials) and implantation plan (pink tacts very close to each other, and c electrodes inserted deeply (pink fiducials) performed with our interface. Moreover, the reason for con- and 14 automatic; power = 90%, p < 0.05, σ = 1 mm) for tacts not being assigned to an electrode was because of bolt the angle error displacement. The sample size has also been heads or contacts not being segmented and because of incor- increased to account for the possible effect of clustering of rectly assigned contacts to bolt heads. However, electrode electrodes within patients, with an assumed ICC = 0.25 and bending did not influence accuracy of contact assignment as average number of 8 electrodes per patient. The mean angle evidenced by our approach performing equivalent between difference observed (paired t test) is small and has a strong the two surgical approaches (placing a stylet close to and far and positive correlation and good reliability. Related to the from TP), where global bending is significantly different. non-inferiority test of the displacement error caused by this angle difference, there is no suggestion that, at the tolerance level of contact length, either method is worse that the other. Validation We were also able to confirm that the anatomy regions at the tip of the electrodes are concordant with those manually The MAE of the centroids validated from the manual iden- identified by two neurosurgeons. Our sample size is above tification of contacts in our study is slightly lower than the the sample size computed (159) that is sufficient for a 95% localisation error of 0.5 mm reported in [2,16], although with confidence interval with width ± 0.1 assuming an estimate of a greater standard deviation. The RMSD reported in our study 0.6 ICC. This is important for post-surgical analysis of SEEG of axial and sagittal planes is similar to the RMS reported in electrodes as knowing the anatomical region each contact is [11]. However, we see a higher error in the coronal plane due located in can aid in identifying the seizure onset zone. to greater CT slice thickness (μ = 0.87 vs. μ = 1.14 mm) and thus a greater RMSD than that reported in their work with deep brain electrodes. We also confirmed the accuracy Bending of the computed contact positions with respect to those from two manual segmentations which varied less that 0.8mm Compared to previous approaches for DBS that fit trajecto- (CI of 95%). We defined equivalence of bolt angles between ries along electrodes using polynomials [10,11], we quantify manual and automatic segmentation as an interval of − 2.29– the amount of bending as well as the displacement at contact 2.29 mm based on the sample size calculation (14 manual positions, permitting to study the reasons of bending within 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 945 the anatomy. The parcellation is used to accurately report the ing a stylet closer to the target point result in lower bending anatomical regions the electrodes have traversed. We were of electrodes. able to estimate local bending by modelling electrodes as Acknowledgements This publication represents in part independent elastic rods and using Darboux vectors to quantify the 3- research commissioned by the Health Innovation Challenge Fund degrees-of-freedom rate of change of the material frames (WT106882), the Wellcome/EPSRC [203145Z/16/Z], and the National orthogonally aligned to the electrode. The large and positive Institute for Health Research University College London Hospitals correlation observed between global bending and the accu- Biomedical Research Centre (NIHR BRC UCLH/UCL High Impact Initiative). We are grateful to the Wolfson Foundation and the Epilepsy mulated displacement of contacts in addition to the medium Society for supporting the Epilepsy Society MRI scanner. The views positive correlation with length of electrode indicate that Dar- expressed in this publication are those of the authors and not necessarily boux vectors are a representative measure of bending. The those of the Wellcome Trust or NIHR. projection of a rigid rod is based on the bolt direction rather than on planned trajectories as in previous studies. This facil- Compliance with ethical standards itates evaluating the displacement of each contact, rather than a displacement due to EP location errors and angle of drilling. Conflict of interest The authors declare that they have no conflict of We confirmed the displacement at the tip of the electrode due interest. to angle difference between bolts automatically and manu- Ethical approval All data were evaluated retrospectively. All studies ally segmented is below a tolerance displacement error, and involving human participants were in accordance with the ethical stan- therefore, we were able to report contact displacement due dards of the institutional and/or national research committee and with to bending with respect to a rigid electrode. the 1964 Helsinki Declaration and its later amendments or comparable ethical standards. Informed consent For this type of study, formal consent is not required. Conclusions and future work Open Access This article is distributed under the terms of the Creative We present a method for automatic segmentation of elec- Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, trodes, including their contacts and bolts, that takes bending and reproduction in any medium, provided you give appropriate credit into account by quantitatively estimating local and global to the original author(s) and the source, provide a link to the Creative bending. We show the importance of accurately detecting Commons license, and indicate if changes were made. the angle of the bolt, since it is one of the main reasons for TP errors, as well as the importance of accurately and auto- matically reporting the region of anatomy the contacts are References located in, since it aids identifying the seizure onset zone. Our approach was validated in 23 data sets comprising two surgi- 1. Arnulfo G, Hirvonen J, Nobili L, Palva S, Palva JM (2015) Phase cal techniques and demonstrated in these cases our method and amplitude correlations in resting-state activity in human stereo- tactical EEG recordings. NeuroImage 15(112):114–127 is robust to bending along the electrode. 2. Arnulfo G, Narizzano M, Cardinale F, Fato MM, Palva JM (2015) Future work is required to guarantee generalisability of Automatic segmentation of deep intracerebral electrodes in com- automatic segmentation of SEEG electrodes by enabling puted tomography scans. BMC Bioinform 16(99):1–12 automatic parameter selection to support data from multi- 3. Cardinale F, Cossu M, Castana L, Casaceli G, Schiariti MP, Mis- erocchi A, Fuschillo D, Moscato A, Caborni C, Arnulfo G, Lo ple centres. We hypothesise that white matter tracks may be Russo G (2013) Stereoelectroencephalography: surgical methodol- one of the factors of electrodes bending, and therefore, we ogy, safety, and stereotactic application accuracy in 500 procedures. envisage using diffusion MRI tractography in combination Neurosurgery 72(3):353–366 with our proposed methods in future studies to understand 4. Cardoso MJ, Modat M, Wolz R, Melbourne A, Cash D, Rueckert D, Ourselin S (2015) Geodesic information flows: spatially-variant the reasons to bending. Understanding the mechanical prop- graphs and their application to segmentation and fusion. IEEE TMI erties of electrodes along with the biomechanical properties 34(9):1976–1988 of the brain tissue as well as simulating instrument–tissue 5. D’Albis T, Haegelen C, Essert C, Fernandez-Vidal S, Lalys F, Jan- interaction will permit greater fidelity to the implantation nin P (2015) PyDBS: an automated image processing workflow for deep brain stimulation surgery. Int J CARS 10(2):117–128 plan resulting in more accurately targeting specific regions 6. Dogdas B, Shattuck DW, Leahy RM (2005) Segmentation of skull and potentially improve clinical outputs including the ability and scalp in 3-D human MRI using mathematical morphology. to reduce the number of implanted electrodes and targeting Hum Brain Mapp 26(4):273–285 riskier areas. We envisage to incorporate our work to an EEG 7. Dorfer C, Minchev G, Czech T, Stefanits H, Feucht M, Pataraia E, Baumgartner C, Kronreif G, Wolfsberger S (2017) A novel minia- analysis pipeline and validate the activity read from SEEG ture robotic device for frameless implantation of depth electrodes contacts with their anatomical location. Parallel clinical work in refractory epilepsy. J Neurosurg 126(5):1622–1628 will look into different types of techniques and their effect on 8. Duncan JS, Sander JW, Sisodiya SM, Walker MC (2006) Adult electrodes bending, i.e. understanding the reasons why push- epilepsy. Lancet 367(9516):1087–1100 123 946 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 9. Duncan JS, Winston GP, Koepp MJ, Ourselin S (2016) Brain 17. Sparks R, Vakharia V, Rodionov R, Vos SB, Diehl B, Wehner imaging in the assessment for epilepsy surgery. Lancet Neurol T, Miserocchi A, McEvoy AW, Duncan JS, Ourselin S (2017) 15(4):420–433 Anatomy-driven multiple trajectory planning (ADMTP) of 10. Husch A, Gemmar P, Lohscheller J, Bernard F, Hertel F (2015) intracranial electrodes for epilepsy surgery. IJCARS 12(8):1245– Assessment of electrode displacement and deformation with 1255 respect to pre-operative planning in deep brain stimulation. Bild- 18. Sparks R, Zombori G, Rodionov R, Nowell M, Vos SB, Zuluaga verarbeitung für die Medizin MA, Diehl B, Wehner T, Miserocchi A, McEvoy AW, Duncan JS, 11. Husch A, Petersen MV, Gemmar P, Goncalves J, Hertel F (2018) Ourselin S (2017) Automated multiple trajectory planning algo- PaCER—A fully automated method for electrode trajectory and rithm for the placement of SEEG electrodes in epilepsy treatment. contact reconstruction in deep brain stimulation. NeuroImage Clin Int J CARS 12(1):123–136 17:80–89 19. Spillmann J, Harders M (2010) Inextensible elastic rods with 12. Kugelstadt T, Schömer E (2016) Position and orientation based torsional friction based on Lagrange multipliers. Comput Anim cosserat rods. In: Eurographics ACM SIGGRAPH symposium on Virtual Worlds 21(6):561–572 computer animation 20. Umetani N, Schmidt R, Stam J (2014) Position-based elastic rods. 13. Lalys F, Haegelen C, D’albis T, Jannin P (2014) Analysis of elec- Eurographics/ACM SIGGRAPH symposium on computer anima- trode deformations in deep brain stimulation surgery. Int J Comput tion pp 1–10 Assist Radiol Surg 9(1):107–117 21. Vakharia VN, Sparks R, O’Keeffe AG, Rodionov R, Miserocchi 14. Meesters S, Ossenblok P, Colon A, Schijns O, Florack L, Boon P, A, McEvoy A, Ourselin S, Duncan J (2017) Accuracy of intracra- Wagner L, Fuster A (2015) Automated identification of intracra- nial electrode placement for stereoencephalography: a systematic nial depth electrodes in computed tomography data. In: IEEE 12th review and meta-analysis. Epilepsia 58(6):921–932 international symposium on biomedical imaging (ISBI) pp 976– 22. van der Loo LE, Schijns OEMG, Hoogland G, Colon AJ, Wagner 979 GL, Dings JTA, Kubben PL (2017) Methodology, outcome, safety 15. Modat M, Cash DM, Daga P, Winston GP, Duncan JS, Ourselin and in vivo accuracy in traditional frame-based stereoelectroen- S (2014) Global image registration using a symmetric block- cephalography. Acta Neurochir 159:1733–46 matching approach. J Med Imaging 1(2):024003 16. Narizzano M, Arnulfo G, Ricci S, Toselli B, Tisdall M, Canessa A, Fato MM, Cardinale F (2017) SEEG assistant: a 3DSlicer extension to support epilepsy surgery. BMC Bioinform 18(124):1–13 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png International Journal of Computer Assisted Radiology and Surgery Springer Journals

Loading next page...
 
/lp/springer_journal/automatic-segmentation-of-stereoelectroencephalography-seeg-electrodes-jg03wOS8RR
Publisher
Springer Journals
Copyright
Copyright © 2018 by The Author(s)
Subject
Medicine & Public Health; Imaging / Radiology; Surgery; Health Informatics; Computer Imaging, Vision, Pattern Recognition and Graphics; Computer Science, general
ISSN
1861-6410
eISSN
1861-6429
D.O.I.
10.1007/s11548-018-1740-8
Publisher site
See Article on Publisher Site

Abstract

Purpose The accurate and automatic localisation of SEEG electrodes is crucial for determining the location of epileptic seizure onset. We propose an algorithm for the automatic segmentation of electrode bolts and contacts that accounts for electrode bending in relation to regional brain anatomy. Methods Co-registered post-implantation CT, pre-implantation MRI, and brain parcellation images are used to create regions of interest to automatically segment bolts and contacts. Contact search strategy is based on the direction of the bolt with distance and angle constraints, in addition to post-processing steps that assign remaining contacts and predict contact position. We measured the accuracy of contact position, bolt angle, and anatomical region at the tip of the electrode in 23 post-SEEG cases comprising two different surgical approaches when placing a guiding stylet close to and far from target point. Local and global bending are computed when modelling electrodes as elastic rods. Results Our approach executed on average in 36.17 s with a sensitivity of 98.81% and a positive predictive value (PPV) of 95.01%. Compared to manual segmentation, the position of contacts had a mean absolute error of 0.38 mm and the mean bolt angle difference of 0.59 resulted in a mean displacement error of 0.68 mm at the tip of the electrode. Anatomical regions at the tip of the electrode were in strong concordance with those selected manually by neurosurgeons, ICC (3, k) = 0.76, with average distance between regions of 0.82 mm when in disagreement. Our approach performed equally in two surgical approaches regardless of the amount of electrode bending. Conclusion We present a method robust to electrode bending that can accurately segment contact positions and bolt orien- tation. The techniques presented in this paper will allow further characterisation of bending within different brain regions. Keywords Epilepsy · SEEG · Automatic segmentation · Bending Introduction Electronic supplementary material The online version of this article (https://doi.org/10.1007/s11548-018-1740-8) contains supplementary Epilepsy is a disease characterised by an enduring predis- material, which is available to authorized users. position to generate epileptic seizures and affects 1% of the population [8]. A third of patients develop chronic refractory B Alejandro Granados focal epilepsy and neurosurgery may provide a cure [9]. alejandro.granados@ucl.ac.uk Brain imaging is fundamental in a typical neurosurgi- Wellcome/EPSRC Centre for Interventional and Surgical cal evaluation for determining the epileptogenic zone (EZ) Sciences, UCL, London, UK National Hospital for Neurology and Neurosurgery, London, The First Affiliated Hospital of Xian Jiaotong University, UK Xian, People’s Republic of China Department of Clinical and Experimental Epilepsy, Institute Vickie and Jack Farber Inst for Neuroscience, Thomas of Neurology, National Hospital for Neurology and Jefferson University, Philadelphia, USA Neurosurgery, London, UK Dementia Research Centre, Department of Neurodegenerative Department of Statistical Science, University College Disease, UCL Institute of Neurology, London, UK London, London, UK 123 936 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 with modalities including structural and functional MRI trodes are modelled fitting a polynomial function and then (e.g. T1/T2-w, FLAIR) and PET [9]. If the EZ is not aligned to a common coordinate system reporting mean identifiable, invasive electroencephalography (EEG) record- deviation that varies from 0.92 to 2.0 mm. However, they ings are performed in the form of stereo-EEG (SEEG) have mostly focused on fitting trajectories using polynomi- or subdural grid insertion. SEEG is a procedure in which als rather than computing the amount of electrode bending multiple electrodes are stereotactically inserted to identify and have not considered the reasons of bending within the the seizure onset zone [21]. Accurate placement of elec- brain anatomy. Although Lalys et al. [13] looked at the rea- trode contacts is important for safety, interpretation of the sons of bending (mainly due to brain shift) by computing recorded electrical signals, and subsequent resection plan- a local and mean curvature index over the entire length of ning [21]. Planning of electrode implantation is crucial for DBS electrodes, the index provides no information about the avoiding blood vessel damage and subsequent intracranial direction of bending. Unlike SEEG procedures, where 8–14 haemorrhage (which occurs in 1–2% of patients), and auto- electrodes are inserted, DBS electrodes are typically inserted matic computer-assisted multiple trajectory planning tools bilaterally and the contacts are very close to the tip. To dis- have been proposed [17,18]. However, intraoperatively, entry criminate between contacts located in white or grey matter, point (EP) accuracy can be affected by misregistration of the Arnulfo et al. [1,16] compute the distance from each contact neuronavigation system, inaccurate alignment, and deflec- to grey–white matter interface. tion during drilling, whereas target point (TP) errors may be caused by the angle at which the electrode passes through skull, deflection of the electrode at the dura or within the Contribution of this paper brain, the rigidity of the electrode, and the depth to which a guiding stylet is inserted [3,21]. Robotic systems have Our main motivation is to automatically segment SEEG con- been introduced to improve EP implantation accuracy [3,7]. tacts and bolts (Ad-Tech Med Instr Corp, USA) relative to However, TP displacement is the main source of error and the anatomy whilst accounting for electrode bending along understanding why and how electrodes bend may help pre- its trajectory at contact positions rather than as a result of TP dict final TP positions during surgical planning and improve displacement. Our algorithm (Fig. 1) allows estimating not EZ localisation [22]. only the position of contacts but also the direction of the bolts Furthermore, it is convenient to have a rapid and reliable inserted into the skull since the angle of the bolt with respect scheme for segmenting contacts, assigning their anatomical to the scalp surface normal is a measure of post-implantation location when interpreting SEEG studies and for guid- accuracy. We quantify local and global bending by means of ing definitive surgical resections. Automatic segmentation electrodes modelled as elastic rods in position-based dynam- approaches have been proposed for SEEG [2,14,16] and deep ics and validate our methods in 23 post-SEEG cases (224 brain stimulation (DBS) [5,10,11] implantation. Arnulfo electrodes, 1843 contacts) comprising two different surgical et al. [2] used post- implantation CT (threshold = 1600) approaches (placing a guiding stylet close to or far from the co-registered with MRI to segment electrodes based on a TP). geometrical-constrained search. They randomly generated different scenarios for 1–15-mm displaced TP in an experi- mental study and reported accuracy of 10% of false negatives Methods (FN) and 7% of false positives (FP) for a maximum dis- placement of 15 mm. However, bending may occur at any Input images point along the electrode’s trajectory. Meesters et al. [14] co-registered the CT (threshold = 500 HU) with MRI and A post-SEEG implantation resampled CT and an MRI T1 extracted guiding screws with a multi-scale filter whilst deter- images are rigidly co-registered using NiftyReg (v1.5.43) mining likely tip locations within a wedge-shape region. [15]. From the MRI image, we obtain the parcellation of However, manual adjustments took between 10 s and several brain anatomy via NiftyWeb (GIF v3.0) (Fig. 2)[4]. minutes, and reported deviations of the tip and their method did not account for electrodes bending. Additionally, these Identification of anatomical masks methods relied on pre-operative plans and were tested only on one electrode type. We use the MRI and the parcellation to create regions of Hubsch et al. [10,11] proposed an automated algorithm interest that are used to identify contacts, bolt heads, and the reconstructing full electrode trajectory whilst accounting for section of the bolt crossing the scalp/skull, which we refer DBS electrode bending from CT scans. A convex hull brain mask is extracted using thresholds, and the largest connected https://github.com/InteractiveComputerGraphics/PositionBasedDyn components are skeletonised [10]. Trajectories of 11 elec- amics. 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 937 Fig. 1 Flow chart of algorithm pipeline Fig. 2 Input images: a post-SEEG implantation CT, b MRI T1, and c parcellation Fig. 3 a Axial, b sagittal and c coronal planes showing computed masks of the brain (cyan), skull and scalp (yellow) together with the result of connected components filters of contacts (red), bolt head (green), and the section of the bolt crossing the skull (blue) as bolt body.First,a BinaryThresholdImageFilter A BinaryThresholdImageFilter is applied to the is applied to the parcellation to create a mask of intracra- MRI to create a mask of the scalp BI with a lower scalp nial space BI , i.e. with a threshold t in the range of threshold equal to t . We use morphological operators brain brain scalp 4 ≤ t ≤ 208. We apply a method similar to Dogdas et al. to combine BI and BI and apply a closing fil- brain brain scalp [6], which we describe herein for completeness. We compute ter with a ball structuring element (radius = 10) to obtain a a skull threshold t from the MRI as the mean of the inten- mask of the head, i.e. BI = (BI ∪ BI )  B , skull head scalp brain 10 sities of the nonzero voxels that are not brain as an empirical and a mask of the skull, i.e. BI = BI ⊕ BI , skull head brain measure to split the low- and high-intensity regions, followed after applying an XOR morphological operator on the result by a scalp threshold t as the mean of the non-brain vox- (Fig. 3). scalp els above the skull threshold (∀I (x , y, z) ≥ t )to MRI MRI skull identify the transition between the head and the background. 123 938 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 Table 1 Geometrical analysis, Geometrical analysis Discriminant analysis μ(σ ), and discriminant analysis Number of Pixels Elongation Roundness Number of pixels Roundness of bolt heads and contacts Bolt head 329.4 (183.5) 2.51 (0.59) 0.63 (0.06) > 100 [0.4, 1.0] Contact 9.7 (6.6) 2.52 (1.27) 1.10 (0.06) [3, 50] Fig. 4 Search strategy given the direction of the bolt and constraints (distance and angle) Segmentation of electrode bolts and contacts (30 ) (Fig. 4), constraints which favour assigning contacts in the direction of the bolt during a first pass. Amask BI is created from a BinaryThreshold post ImageFilter applied to the post-op CT with lower thresh- Automatic segmentation of electrodes old t = (0.52) ∗ max (I (x , y, z)). BI is used to CT CT post identify full bolts (BI ) with at least a minimum of 200 bolt The main steps of our algorithm include: pixels. Three subsections are identified: the head of the bolt which is outside the patient’s head (BI ∩¬BI ), the 1. Initialisation All segmented contacts are initially labelled bolt head as ’available’ and stored in a pool. Given a bolt head body (BI ∩ BI ), i.e. section crossing the skull, and the bolt skull position (x ), the closest bolt bodies (8 ≤ x − x ≤ tip (BI ∩ BI ). Lastly, contacts are identified within h h b bolt brain 25 mm) and the closest contact ( x − x ≤ 50 mm) the brain whilst excluding bolt tips ((BI ∩ BI ) ⊕ h c post brain are identified in order to narrow the search down to only BI ). We applied a ConnectedComponentImage boltTip Filter to the masks and a LabelImageToShape those relevant. 2. Contact search strategy For each bolt head, the contact LabelMapFilter to the blobs to get their centroids and geometrical properties before conducting geometrical anal- search strategy is executed initially with the closest bolt body (1st pass search) and subsequently with alternative ysis to identify discriminants of segmentation (Table 1). We detected contacts with blobs that were within a range of num- bolt bodies if no contacts have been assigned. Although rare, bolt bodies may not be segmented and a direction of ber of pixels ([3, 50]) and bolt heads with blobs that had a minimum number of pixels (≥ 100) and were within a range search cannot be computed. Therefore, the contact search strategy is called again with the closest contact position of roundness values ([0.4, 1.0]). rather than a bolt body position. 3. Project remaining contacts in pool For electrodes con- Contact search strategy taining at least one contact, we compute the minimum distance between an available contact in the pool and Given a bolt head (x ) and its closest bolt body (x ) posi- a line formed by the positions of the bolt head and the h b x −x b h tions, we compute the direction of search (nˆ = ) electrode tip. The contact is assigned to the electrode if x −x b h and iteratively compute a number of points x given a maxi- and only if its distance to the closest point x to the line p p mum electrode length (90 mm) and a step size (1 mm) in the (tangent to the line) is below a constraint (5 mm) and x direction nˆ . An available contact x is assigned to the elec- remains along the line or in a position of the line 20% 0 c trode if and only if it is located below a distance constraint extended from the tip, i.e. within an interpolation range from x (5 mm) and the angle between the previous direction of [0.0, 1.2] to project contacts that are further from the nˆ and the current direction nˆ is below an angle constraint currently identified tip of the electrode. 0 c 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 939 Fig. 5 Modelling of electrodes as elastic rods. Bolt head (green) and body (blue) with contacts (red) modelled as point particles and ghost particles (cyan) created orthogonally along the electrode with material frames located between contacts 4. Predict contacts in the bolt region For a given electrode, electrodes (N = 109 contacts), (b) manually identify the we compute the most common segment along the elec- tip and head of the bolt of a random subset of electrodes trode based on the distances between subsequent contacts (N = 95 bolts), (c) confirm the correct number and location rounded to the closest integer. Based on electrode spec- of contacts and electrodes (N = 23 cases), and (d) identify ification, we infer the type of electrode depending on the TP anatomical region (N = 222 electrodes). the order of the segments and specify contact spacing. We then compute the direction from the last contact x towards the bolt head x and create new contacts up to 21 Results mm before the bolt head position to segment only those contacts closer to the skull. Interface Bending estimation We implemented our algorithms in C++ using MITK and ITK as well as a GUI in Qt to allow clinicians to adjust the To quantify electrode bending, electrodes are modelled as automatic segmentation if needed (Fig. 6). On average, our elastic rods using the Cosserat model proposed by [19] method executed in 36.17 s (N = 23, σ = 15.7), faster than and then implemented by [12] in position-based dynam- manual segmentation. ics. Electrode contact positions are represented as linked particles with ghost particles located orthogonally half-way between contact pairs (Fig. 5). A material frame is created Performance between contacts with a unit vector (d = X − X ) 3 c c n−1 n aligned tangentially to its centreline followed by two addi- Of a total of 224 electrodes (1843 contacts), 29 contacts were tional orthonormal vectors, (d = d × (X − X )) and 2 3 c c n−1 segmented but not assigned to any electrode due to: (a) three ˆ ˆ (d = d × d ) chosen to lie in the principal direction of the 1 2 3 bolt heads that were not automatically segmented (17), (b) cross section. We compute the rate of change of two consec- no segmented contacts close to them (5), and (c) due to one utive frames, namely a Darboux vector , to describe local incorrectly assigned contact to a bolt head (7). On average, bending at the contact points [12,20]. Along the electrode, Tr P Tr P the sensitivity ( ∗ 100) and PPV ( ∗ 100) of Tr P+FN Tr P+FP values are then accumulated to quantify global bending. We our approach was μ = 98.81%; σ = 2.04 (false-negative then use the parcellation to report the region at which each rate of μ = 0.124; σ = 0.02) and μ = 95.01%; σ = 6.73 contact is located and report all those regions that the elec- (false- positive rate of μ = 0.059; σ = 0.09), respectively trode passes through. Lastly, contact displacement and depth (Fig. 7 bottom), finding no statistical significant difference are estimated with respect to a rigid electrode with position between data sets of the two surgical approaches, i.e. placing of contacts projected along the direction from the bolt head to a stylet far from or close to target point. To illustrate our the last contact (X ) at distances subject to electrode speci- results, Fig. 8 shows two cases correctly identified (a, b) fication. along two worst cases (c, d). Validation We asked two neurosurgeons and one clinical scientist to http://mitk.org. (a) manually segment the contacts of a random subset of http://itk.org/. 123 940 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 Fig. 6 Automatic segmentation of electrodes interface and GUI for manual adjustment Validation Computed contact positions, bolt angles, and regions of anatomy are compared with those manually segmented in a subset of cases (Table 2). – Contact position Compared to the manual segmentation done by a clinical scientist (M1) and a neurosurgeon (M2), we found that the contact location of our auto- matic segmentation approach had a mean absolute error (MAE) of 0.38 and 0.40 mm, respectively, and a root- mean-square deviation (RMSD) of 0.45. The distance of contact positions between both manual segmentations was on average μ = 0.37 mm (σ = 0.22). We found no statistical difference when comparing the distances from automatically computed contact position to those positions obtained via manual segmentation (paired dif- ferences: μ = 0.036, σ = 0.21). – Bolt angle We found that the angle of bolts between automatic and manual segmentation (Fig. 9) by M1 and ◦ ◦ M2 differed on average by 0.59 and 0.22 , respec- tively, with pair samples strongly and positively corre- lated (Pearson correlation) and with strong reliability (Cronbach’s alpha). We study the displacement error d = sin(θ ) ∗ l at the tip of the electrode error error e caused by this angle difference θ and the length error of the electrode l within the brain and define a maxi- mum tolerance value T = 2.29 mm related to contact length. On average, d at the tip of a rigid elec- error trode caused by the angle difference is μ = 0.68 mm Fig. 7 Top: Number of contacts initially segmented and assigned to for M1 and μ = 0.72 mm for M2. We found 3 out- electrodes via a bolt head and bolt bodies association (step 2), pro- jected (step 3), predicted (step 4), and left unassigned and available liers above T for M1 ([2.46, 4.79] mm) and 8 outliers in pool. Bottom: Number of contacts correctly segmented (TrP—true for M2 ([2.37, 5.48] mm). Given T , a non-inferiority positives), wrongly segmented (FP) and missed (FN) in 23 data sets test indicates that 0.68 mm is an estimate of d with error (order by number of electrodes) 95% of CI (0.431–0.926) after accounting for clustering 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 941 Fig. 8 Examples: a segmented bolt head and contacts of electrodes overlaying CT; b contacts predicted at the skull and scalp level; c 22 FP (red marks along the skull); and d our worst case with 3 contacts not segmented due to crossings and 4 FNs using a patient-level random effect. Figure 10 shows of electrodes (μ = 0.49; σ = 0.34) compared to the bend- an example with electrodes automatically segmented ing observed in electrodes that had a stylet placed close to TP and their corresponding rigid electrodes (lighter colours) (μ = 0.31; σ = 0.18), a difference which was statistically computed using the direction of bolts automatically significant, t (222) = 5.36, p < 0.01. segmented. – Regions of anatomy We also ran a intra-class correlation Generalisability and robustness two-way mixed effects model with average measures and found a strong agreement when identifying the anatomi- Three SEEG post-resection cases using SEEG DEPTH elec- cal region of the brain at the tip of the electrode between trodes (PMT Corp., USA) were obtained from the Vickie our algorithm and that done manually by two neuro- and Jack Farber Institute for Neuroscience (Thomas Jeffer- surgeons, ICC (3, k) = 0.76, p < 0.001. When in son University) to assess generalisability and robustness of disagreement, the average distance between regions was our algorithm. We observed a reduced average performance 0.82 mm (σ = 0.78), a distance below contact size. (sensitivity = 69.7% and PPV = 82.6%) due to the following Furthermore, it is estimated that electrode contacts elec- factors (Fig. 11a–c): (a) smaller bolt heads (our parameter of trically sample regions of grey matter within a 3 mm minimum number of pixels of bolt heads could be adjusted), radius. Any discrepancy in identified anatomical regions (b) contacts being very close to each other and merged as sin- below this is therefore not clinically significant. gle blobs (addressed by adopting optimal oblique resampling used for DBS electrodes [11]), and (c) electrodes inserted Bending estimation deeply (our parameter of maximum electrode length could be adjusted to account for this). Despite this, our algorithm In order to study whether Darboux vectors are a represen- was agnostic of electrode types and implantation planning tative measure of bending, we look into the relationship and was robust in post-resection data sets. (Pearson correlation) between global bending and the fol- We randomly chose 3 of our data sets to test the method lowing variables: accumulated displacement of contacts (r = proposed in [2,16] and implemented in SEEGA (Slicer 0.532, p < 0.001), length of electrode inside the brain tissue v4.6.2). We configured electrode types based on electrode (r = 0.373, p < 0.001), amount of white matter traversed specification, used the implantation plan (EP and TP) as fidu- by the electrode (r = 0.257, p < 0.001), and bolt angle cials and imported the CT image. We modified SEEGA to (r = 0.189, p = 0.045). Of the two surgical approaches, use the same threshold that we computed in our algorithm placing a stylet far from TP resulted in larger global bending for consistency and because the default threshold computed 123 942 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 Table 2 Validation between manual and automatic segmentation Measure Manual vs automatic Statistical test Result Contact position MAE (μ, σ , IQR) M1: 0.38 mm, 0.24, 0.22 Paired t test t (106) = − 1.756, p = 0.82 N = 109/1843 M2: 0.40 mm; 0.22, 0.26 Pearson correlation r = 0.454, p < 0.001 Cronbach’s alpha 0.615 MAE (x , y, z M1: (0.14, 0.15, 0.27) mm components) M2: (0.17, 0.15, 0.26) mm RMSD M1: 0.45 M2: 0.45 Bolt angle mean angle difference M1: 0.59 (1.27) Paired t test t (94) = − 4.54, p < 0.001 M1: N=95/224 M2: 0.22 (1.53) t (112) = 1.533, p = 0.128 Pearson correlation r = 0.991, p < 0.001 r = 0.985, p < 0.001 Cronbach’s alpha 0.995 0.992 M2: N = 113/224 displacement error at M1: 0.68 mm, 0.81, 0.83 Non-inferiority test CI = (0.431, 0.926) first contact due to tolerance = 2.29 mm M2: 0.72 mm, 0.84, 0.81 angle difference (μ, σ,IQR) Regions of anatomy region of anatomy at Intra-class correlation 0.76, p < 0.001 N = 222/222 first contact distance between 0.82 (0.78) mm regions when in disagreement Fig. 9 Bolt angles. a Bolt from post-CT image and b manual iden- fication (rigid electrode shown in blue) of bolt direction of an outlier tification of the direction along bolts by a clinical scientist. c, d case with angle difference of θ = 5.73 and displacement error at error (Inconspicuous) comparison of manual (pink) and automatic identi- the tip of d = 4.79 mm error by SEEGA resulted in segmentations errors. We observed noticed that most of the contact positions were not accurate an average sensitivity of 82.9% and PPV of 65.3% (97.3% (Fig. 11d), we only report performance on the number of con- and 98.2%, respectively, using our algorithm). Whilst we tacts that were incorrectly or not segmented. Similarly to our 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 943 Fig. 10 Anatomical regions traversal (top) and contact displacement (bottom) of automatically segmented electrodes with respect to a rigid electrode computed based on bolt direction approach tested with data from other centres, SEEGA might We found that FP were located in the inner surface of the perform better after fine-tuning parameters. Both algorithms skull and were caused by pixel size inaccuracies of CSF have parameters the user must select to guarantee optimal regions overlapping with bone structure. The performance performance. of our algorithm is similar to previous approaches although [2] only considers displacements at the tip of the electrode with no details of displacement of other contacts along the Discussion electrode and [14] uses a very small sample size assuming rigid electrodes. Automatic segmentation Compared to the search strategy by Arnulfo et al. [2], our algorithm uses a higher angle constraint (30 rather than The automatic segmentation of bolts is typically overlooked 10 ) because we use instead the bolt direction to search for in the literature and could be used to report accuracy errors contacts rather than a direction from previously segmented caused by differences in angle with respect to planning. contacts. It is also clinically relevant to accurately segment We use bold direction to search for contacts with neither the position of the contacts closer to the skull to ensure grey prior information of electrode type nor implantation plan- matter at the cortical entry is adequately sampled, but these ning. Compared to previous work [2,14], we use a factor might be difficult to segment in the bolt region. Therefore, of maximum intensity from CT images rather than a con- contacts are predicted in this region after inferring electrode stant. Regions of interests were used to segment position type. Further conditions would need to be included in this step of contacts and bolts based on geometrical properties with to support more electrodes from different manufacturers. To their centroids equivalent to the signal peaks found in [11] cope with electrodes crossing, we initially used geometri- and more generally in the literature. The choice of intensity cal features to identify large blobs that relate to more than threshold and constraints favour few incorrectly segmented one contact for splitting. However, the resulting position of contacts (FP) over missing contacts (FN), since these can contacts was not good enough to make the method fully auto- be easily discarded by surgeons during manual adjustment. matic, so we rely on manual adjustments that can be quickly 123 944 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 Fig. 11 Generalisability and robustness tests. Our proposed algorithm electrode with insertion depth of 110 mm); Our data in SEEGA: d seg- using data from a different centre: a with smaller bolt heads, b con- mented contact positions (green fiducials) and implantation plan (pink tacts very close to each other, and c electrodes inserted deeply (pink fiducials) performed with our interface. Moreover, the reason for con- and 14 automatic; power = 90%, p < 0.05, σ = 1 mm) for tacts not being assigned to an electrode was because of bolt the angle error displacement. The sample size has also been heads or contacts not being segmented and because of incor- increased to account for the possible effect of clustering of rectly assigned contacts to bolt heads. However, electrode electrodes within patients, with an assumed ICC = 0.25 and bending did not influence accuracy of contact assignment as average number of 8 electrodes per patient. The mean angle evidenced by our approach performing equivalent between difference observed (paired t test) is small and has a strong the two surgical approaches (placing a stylet close to and far and positive correlation and good reliability. Related to the from TP), where global bending is significantly different. non-inferiority test of the displacement error caused by this angle difference, there is no suggestion that, at the tolerance level of contact length, either method is worse that the other. Validation We were also able to confirm that the anatomy regions at the tip of the electrodes are concordant with those manually The MAE of the centroids validated from the manual iden- identified by two neurosurgeons. Our sample size is above tification of contacts in our study is slightly lower than the the sample size computed (159) that is sufficient for a 95% localisation error of 0.5 mm reported in [2,16], although with confidence interval with width ± 0.1 assuming an estimate of a greater standard deviation. The RMSD reported in our study 0.6 ICC. This is important for post-surgical analysis of SEEG of axial and sagittal planes is similar to the RMS reported in electrodes as knowing the anatomical region each contact is [11]. However, we see a higher error in the coronal plane due located in can aid in identifying the seizure onset zone. to greater CT slice thickness (μ = 0.87 vs. μ = 1.14 mm) and thus a greater RMSD than that reported in their work with deep brain electrodes. We also confirmed the accuracy Bending of the computed contact positions with respect to those from two manual segmentations which varied less that 0.8mm Compared to previous approaches for DBS that fit trajecto- (CI of 95%). We defined equivalence of bolt angles between ries along electrodes using polynomials [10,11], we quantify manual and automatic segmentation as an interval of − 2.29– the amount of bending as well as the displacement at contact 2.29 mm based on the sample size calculation (14 manual positions, permitting to study the reasons of bending within 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 945 the anatomy. The parcellation is used to accurately report the ing a stylet closer to the target point result in lower bending anatomical regions the electrodes have traversed. We were of electrodes. able to estimate local bending by modelling electrodes as Acknowledgements This publication represents in part independent elastic rods and using Darboux vectors to quantify the 3- research commissioned by the Health Innovation Challenge Fund degrees-of-freedom rate of change of the material frames (WT106882), the Wellcome/EPSRC [203145Z/16/Z], and the National orthogonally aligned to the electrode. The large and positive Institute for Health Research University College London Hospitals correlation observed between global bending and the accu- Biomedical Research Centre (NIHR BRC UCLH/UCL High Impact Initiative). We are grateful to the Wolfson Foundation and the Epilepsy mulated displacement of contacts in addition to the medium Society for supporting the Epilepsy Society MRI scanner. The views positive correlation with length of electrode indicate that Dar- expressed in this publication are those of the authors and not necessarily boux vectors are a representative measure of bending. The those of the Wellcome Trust or NIHR. projection of a rigid rod is based on the bolt direction rather than on planned trajectories as in previous studies. This facil- Compliance with ethical standards itates evaluating the displacement of each contact, rather than a displacement due to EP location errors and angle of drilling. Conflict of interest The authors declare that they have no conflict of We confirmed the displacement at the tip of the electrode due interest. to angle difference between bolts automatically and manu- Ethical approval All data were evaluated retrospectively. All studies ally segmented is below a tolerance displacement error, and involving human participants were in accordance with the ethical stan- therefore, we were able to report contact displacement due dards of the institutional and/or national research committee and with to bending with respect to a rigid electrode. the 1964 Helsinki Declaration and its later amendments or comparable ethical standards. Informed consent For this type of study, formal consent is not required. Conclusions and future work Open Access This article is distributed under the terms of the Creative We present a method for automatic segmentation of elec- Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, trodes, including their contacts and bolts, that takes bending and reproduction in any medium, provided you give appropriate credit into account by quantitatively estimating local and global to the original author(s) and the source, provide a link to the Creative bending. We show the importance of accurately detecting Commons license, and indicate if changes were made. the angle of the bolt, since it is one of the main reasons for TP errors, as well as the importance of accurately and auto- matically reporting the region of anatomy the contacts are References located in, since it aids identifying the seizure onset zone. Our approach was validated in 23 data sets comprising two surgi- 1. Arnulfo G, Hirvonen J, Nobili L, Palva S, Palva JM (2015) Phase cal techniques and demonstrated in these cases our method and amplitude correlations in resting-state activity in human stereo- tactical EEG recordings. NeuroImage 15(112):114–127 is robust to bending along the electrode. 2. Arnulfo G, Narizzano M, Cardinale F, Fato MM, Palva JM (2015) Future work is required to guarantee generalisability of Automatic segmentation of deep intracerebral electrodes in com- automatic segmentation of SEEG electrodes by enabling puted tomography scans. BMC Bioinform 16(99):1–12 automatic parameter selection to support data from multi- 3. Cardinale F, Cossu M, Castana L, Casaceli G, Schiariti MP, Mis- erocchi A, Fuschillo D, Moscato A, Caborni C, Arnulfo G, Lo ple centres. We hypothesise that white matter tracks may be Russo G (2013) Stereoelectroencephalography: surgical methodol- one of the factors of electrodes bending, and therefore, we ogy, safety, and stereotactic application accuracy in 500 procedures. envisage using diffusion MRI tractography in combination Neurosurgery 72(3):353–366 with our proposed methods in future studies to understand 4. Cardoso MJ, Modat M, Wolz R, Melbourne A, Cash D, Rueckert D, Ourselin S (2015) Geodesic information flows: spatially-variant the reasons to bending. Understanding the mechanical prop- graphs and their application to segmentation and fusion. IEEE TMI erties of electrodes along with the biomechanical properties 34(9):1976–1988 of the brain tissue as well as simulating instrument–tissue 5. D’Albis T, Haegelen C, Essert C, Fernandez-Vidal S, Lalys F, Jan- interaction will permit greater fidelity to the implantation nin P (2015) PyDBS: an automated image processing workflow for deep brain stimulation surgery. Int J CARS 10(2):117–128 plan resulting in more accurately targeting specific regions 6. Dogdas B, Shattuck DW, Leahy RM (2005) Segmentation of skull and potentially improve clinical outputs including the ability and scalp in 3-D human MRI using mathematical morphology. to reduce the number of implanted electrodes and targeting Hum Brain Mapp 26(4):273–285 riskier areas. We envisage to incorporate our work to an EEG 7. Dorfer C, Minchev G, Czech T, Stefanits H, Feucht M, Pataraia E, Baumgartner C, Kronreif G, Wolfsberger S (2017) A novel minia- analysis pipeline and validate the activity read from SEEG ture robotic device for frameless implantation of depth electrodes contacts with their anatomical location. Parallel clinical work in refractory epilepsy. J Neurosurg 126(5):1622–1628 will look into different types of techniques and their effect on 8. Duncan JS, Sander JW, Sisodiya SM, Walker MC (2006) Adult electrodes bending, i.e. understanding the reasons why push- epilepsy. Lancet 367(9516):1087–1100 123 946 International Journal of Computer Assisted Radiology and Surgery (2018) 13:935–946 9. Duncan JS, Winston GP, Koepp MJ, Ourselin S (2016) Brain 17. Sparks R, Vakharia V, Rodionov R, Vos SB, Diehl B, Wehner imaging in the assessment for epilepsy surgery. Lancet Neurol T, Miserocchi A, McEvoy AW, Duncan JS, Ourselin S (2017) 15(4):420–433 Anatomy-driven multiple trajectory planning (ADMTP) of 10. Husch A, Gemmar P, Lohscheller J, Bernard F, Hertel F (2015) intracranial electrodes for epilepsy surgery. IJCARS 12(8):1245– Assessment of electrode displacement and deformation with 1255 respect to pre-operative planning in deep brain stimulation. Bild- 18. Sparks R, Zombori G, Rodionov R, Nowell M, Vos SB, Zuluaga verarbeitung für die Medizin MA, Diehl B, Wehner T, Miserocchi A, McEvoy AW, Duncan JS, 11. Husch A, Petersen MV, Gemmar P, Goncalves J, Hertel F (2018) Ourselin S (2017) Automated multiple trajectory planning algo- PaCER—A fully automated method for electrode trajectory and rithm for the placement of SEEG electrodes in epilepsy treatment. contact reconstruction in deep brain stimulation. NeuroImage Clin Int J CARS 12(1):123–136 17:80–89 19. Spillmann J, Harders M (2010) Inextensible elastic rods with 12. Kugelstadt T, Schömer E (2016) Position and orientation based torsional friction based on Lagrange multipliers. Comput Anim cosserat rods. In: Eurographics ACM SIGGRAPH symposium on Virtual Worlds 21(6):561–572 computer animation 20. Umetani N, Schmidt R, Stam J (2014) Position-based elastic rods. 13. Lalys F, Haegelen C, D’albis T, Jannin P (2014) Analysis of elec- Eurographics/ACM SIGGRAPH symposium on computer anima- trode deformations in deep brain stimulation surgery. Int J Comput tion pp 1–10 Assist Radiol Surg 9(1):107–117 21. Vakharia VN, Sparks R, O’Keeffe AG, Rodionov R, Miserocchi 14. Meesters S, Ossenblok P, Colon A, Schijns O, Florack L, Boon P, A, McEvoy A, Ourselin S, Duncan J (2017) Accuracy of intracra- Wagner L, Fuster A (2015) Automated identification of intracra- nial electrode placement for stereoencephalography: a systematic nial depth electrodes in computed tomography data. In: IEEE 12th review and meta-analysis. Epilepsia 58(6):921–932 international symposium on biomedical imaging (ISBI) pp 976– 22. van der Loo LE, Schijns OEMG, Hoogland G, Colon AJ, Wagner 979 GL, Dings JTA, Kubben PL (2017) Methodology, outcome, safety 15. Modat M, Cash DM, Daga P, Winston GP, Duncan JS, Ourselin and in vivo accuracy in traditional frame-based stereoelectroen- S (2014) Global image registration using a symmetric block- cephalography. Acta Neurochir 159:1733–46 matching approach. J Med Imaging 1(2):024003 16. Narizzano M, Arnulfo G, Ricci S, Toselli B, Tisdall M, Canessa A, Fato MM, Cardinale F (2017) SEEG assistant: a 3DSlicer extension to support epilepsy surgery. BMC Bioinform 18(124):1–13

Journal

International Journal of Computer Assisted Radiology and SurgerySpringer Journals

Published: May 7, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off