Home

This is a basic tutorial for the functionality in ielu.

Overview

ielu stands for interactive electrode localization utility. Previously, the project was called gselu stands for grids and strips electrode localization utility, for which the git repository is named.

The purpose of the program is to semi-automatically sort and extract the electrode locations for surface grids and strips, as well as for depth leads.

However, a lot of the interesting algorithms that exist to determine the structures hit from the trajectory of the depth electrodes, are not implemented at all in gselu, reflecting its emphasis primarily on grid-type electrodes.

Requirements and installation

When using ielu, I strongly recommend to use a scientific python distribution such as anaconda or canopy. We primarily use anaconda.

If you are rolling your own distribution here are the dependencies you need:

• traitsui
• chaco
• mayavi / tvtk

And here are the dependencies included in the installation package:

• pysurfer
• nibabel
• mne-python
• pymcubes

These are all open-source libraries. You can easily install them with tools such as pip and easy-install.

It may be helpful to set the environment variables ETS_TOOLKIT='qt4', as well as QT_API which can be either 'pyqt' or 'pyside'. By default, canopy uses pyside, but this occasionally causes dependency problems.

Martinos notes

• The martinos center's python installation does not have pymcubes. It is recommended to roll your own python distribution.

Basic usage

Gselu basically takes two inputs: a CT image (in an image format nibabel can read such as .nii.gz or .mgz), and a freesurfer reconstruction (Note that gselu needs write permission to the freesurfer reconstruction as it will create some surface transformation files. It will not overwrite anything).

Shown below is the main GUI. The CT image is specified in the box labeled "Ct scan" and the freesurfer reconstruction is specified in the boxes labeled "Subjects dir" and "Subject." The SUBJECTS_DIR and SUBJECT environment variables can alternately be set in the underlying environment and these boxes will be populated automatically.

The main pipeline in gselu is accessed by clicking "run pipeline". What this does is:

• Automatically extract the candidate electrode locations from the CT image.
  • This is done by applying a threshold (which is set to a default intensity value but can and should be change for each CT in the parameters menu) and then clustering them by applying a simple spherical binary erosion kernel. The center of mass for each cluster marks each candidate electrode location.
• Automatically register the CT and structural MR images
  • This uses freesurfer's mri_robust_register which has an undocumented option to use minimize normalized mutual information as its global cost function. NOTE - This does not always work perfectly, see the section on registration errors below.
• Use the registration matrix generated above to register the brain mask to the CT image and discard any candidate electrodes outside of the brain.
• Automatically sort the electrodes into candidate grids and strips by using a priori information about the grid geometry.
  • The algorithm to sort the electrodes into geometrically plausible grids incorporates information about the grids' local geometry -- that is, it expects an rectangular NxM grid with a fixed distance f between electrodes. (The algorithm deduces the value of f, but the user provides the values of N and M The algorithm iterates through the candidate electrode locations and greedily adds electrodes to the grid whenever they satisfy the distance and angle conditions, with some tolerance parameters.
  • The grid geometry must be specified in the "Electrode geometry" boxes on the GUI, as shown above. Clicking on the arrows will allow you to add more or fewer grids to this list.
  • The tolerance parameters can be adjusted in the options menu.
  • The algorithm is greedy, and removes points assigned to one grid from further use. So bigger grids should be listed first.
  • This algorithm only applies to subdural grids of size NxM, where N, M >= 2. Depth electrodes or Nx1 subdural strips should not be indicated in the input geometry, they can be assigned manually afterwards.

Following the completion of the base pipeline, some subsequent steps require some user feedback, or can be optionally performed.

• Visualize the sorted grids and strips using mayavi and pysurfer.
  • Visualizations of the electrodes in the space of the surface, in snapped-to-surface coordinates, and in the space of the CT image are shown. The geometric relationships can sometimes be more easily seen in the space of the CT image.
• Edit the sorted grids and strips.
  • In case the algorithms for automatically sorting and extracting the grids/strips did not work perfectly, the user can interactively adjust the solution.
• Automatically label the sorted grids and strips.
  • The sorting algorithm is re-run to associate the electrodes with channel names as found in the raw data. The user should correct rotations or reflections of the correct orientation, as depicted on the surgical schematics.
• Use freesurfer's conventions to translate the sorted electrodes to the surface space, then snap these electrodes to plausible surface coordinates.
  • This uses an iterative algorithm to snap electrodes to surface vertices. It is described in Dykstra et al. 2012. This step is only relevant to subdural surface electrodes, rather than depth leads.

An example of the completed pipeline is shown below:

Note that the version of the GUI shown above is an older version of the software, and some of the buttons and menus have changed.

As you can see, for this subject the procedure works pretty well. However, there are some minor errors. While the software automatically tries to make corrections for these errors, the main purpose of our program is to find a reasonably good solution quickly, and then provide a powerful interactive GUI to allow the user to quickly identify and correct them.

Assigning electrodes to depth leads

Depth electrodes which have Nx1 geometry can't be detected by the sorting and labeling algorithm.

There are some algorithms that exist for identifying the trajectory of depth leads based on the extracranial scattering artifact caused by the leads. In the long term I am considering to apply some of these types of algorithms into ielu eventually, but it is difficult and not the main focus of the software.

Instead depth electrodes should be assigned manually. To do so, the user should click on "add new grid". This tells ielu that there is another lead/grid/strip which has not been previously defined. The software refers to all of these groupings as "grids."

A new grid using a distinct color will appear in the grid menu, shown as follows:

The user selects a particular grid from the drop down menu, and can then click on the electrodes in the 3D view belonging to this lead/grid/strip. They will automatically change color and be added to the program's internal representation of the grid.

This functionality is best demonstrated with a short video, although note that the video was made with an older version of the software and all the buttons in different places:

https://www.youtube.com/watch?v=E2Ol9snQKEk

Note - The first several colors chosen when the user clicks "Add new grid" are hard-coded and provide high contrast between one another. At some point, continuing to add new grids will cause the new grids to have a randomly selected color. It's fairly hard to make this color be dynamically adjustable although I might consider doing so in the future. If you encounter a color that is not a good color for any reason, there is no penalty to ignoring that grid, assigning no electrodes to it, and simply creating and using a different one).

Editing grids

The sorting and labeling algorithm assigns electrodes to best-fit grids. In some minority of cases, these assignments are wrong and the user has to correct them. The process of making these corrections is very similar to the process described above, of identifying depth leads.

First, navigate to the grid in question from the drop-down grid menu. This causes you to assign or remove electrodes from that particular grid. Then click on any electrodes whose grid assignments need to be changed.

Creating new electrodes

In some cases the grid assignment process might require you to assign some electrodes to a position where there is currently a blank space.

There are two ways to create new electrodes at a point where the software doesn't currently know there should be an electrode. These are manual identification from image and linear interpolation. Manually adding the electrodes tends to be a little bit easier.

Manually adding electrodes from CT

From the main GUI, select the menu item Tools -> Examine CT. After a moment, a window displaying the CT image in a 3-slice view should appear.

• Left clicking will move the cursor around the display.
• If there is a pin (as shown below), right clicking will move it to a new position. You don't need to use a pin in order to place a new electrode. A pin only currently appears if from the electrode menu, you do Operations -> Manually modify electrode position. The pin tracks the electrode selected from the electrode menu.

When you are satisfied that you have the cursor or pin in to the right place, you can make one of several selections:

• Track cursor : Adds a white tracking marker to the display, to assist with identifying the current cursor location. Clicking 'track cursor' a second time will remove the tracking marker.
• Add new electrode here : Adds an electrode at the cursor's current location.
• Move electrode : This option moves the electrode represented by the current pin to the pin's new location. If snapping has been performed, clicking this button will undo the snapping.
• Move for postprocessing : This option moves the electrode represented by the current pin to the pin's new location, after snapping. That is, the program will represent the electrode's new location as its "after snapping" location, rather than its normal location. Therefore this option exists just to tweak electrode locations for pretty-picture purposes. ROIs calculated to be in proximity to this electrode will not be altered (no matter how far the electrode has moved).

######Linear interpolation The steps for linear interpolation are as follows:

1. Create a new electrode using the menu item "Create blank electrode". This electrode has no location / grid information, which we will fill in using the rest of the grid's geometry.

2. Fill out the geometry for every electrode manually, as well as the geometry of the electrode that is missing. That is to say, in order to position an electrode using linear interpolation, we need to know which electrodes are immediately surrounding it. So for instance, if we were to place the electrode at point (4,3) we would need to know either where are (4,2) and (4,4), or it would suffice to know (3,3) and (5,3). Or even, it would suffice to know (2,3) and (5,3). In short, it's ok if the geometry information is underspecified, as long as it is logically sufficient to work out a linear interpolation.

   1. Note - it's typically easiest to calculate the geometry on whatever grid the pipeline returns as the best fit, before making corrections, and then make manual adjustments if the geometry returned has some errors.
   2. Note - the geometry uses a CSVListEditor for filling in grids, and it is a little sensitive. Including brackets of any kind ([)] will result in the GUI rejecting your input. It expects a comma-separated list of numbers, with no brackets.

3. When the geometry coordinates are fully specified, select the new electrode and click on "Linearly interpolate". This will cause the electrode to be positioned and the window to reload.

Note: Clicking the linearly interpolate menu item only interpolates one electrode at a time; the currently selected electrode. If you close the electrode window, any "blank" electrodes that haven't been interpolated yet, will be thrown away.

####Electrode window Click on the "examine electrodes" button to open the electrode window. You can do quite a few things in this window.

This window shows the electrodes that have been assigned, whether automatically or manually, to an individual grid -- the grid currently selected in the drop down grid menu. Everything you can do from here, only applies to a single grid. Note, that the electrodes are locked at the time the menu is opened; no more electrodes (except for linear interpolations) can be included. To make changes, close the window, make the changes, and then open it again.

Also, to ensure that bad buggy behavior doesn't happen, you can only open one instance of the electrode window at a time. To do some operations on a different grid/lead/strip, requires to finish what you are doing, close the electrode window, switch to the other grid in the grid menu, and then open the electrode window up again.

#####Labeling - Manual or automatic

To associate the electrodes from a given grid with the channel names provided in the raw data file, first open the electrode window (By clicking the "Examine electrodes" button)

One way to assign names to the electrodes is to manually enter in the names, as is done in the video above, in an earlier version of the software.

However, since that video was created we have additionally developed functionality which automatically tries to determine the names, without user input.

If you select the "Label automatically" button in the electrode window or the menu action of the same name (these do the same thing), the program will automatically try to determine the geometry configuration from the constrained set of points now included in the grid. This procedure should be done after removing any incorrectly sorted electrodes and manually replacing them with the correctly sorted electrodes.

There is another menu item which works by fitting the points onto a plane and optimizing a fast constraint problem to determine their positions in the grid. It doesn't work that well but we include it in the Labeling menu as an option in case the default labeling algorithm fails. Unlike the default labeling algorithm, the plane fitting version is guaranteed to terminate, but might end up with some errors in the result.

The "naming convention" parameter has several possible values

• grid_serial - The electrode in position (7,8) would be called, for instance, G56 since (7-1)*8 + 8 = 56. This is probably the most common naming convention.
• grid_concatenate - This naming convention is for small grids. The electrode in position (7,8) would be called, for instance, G78
• line - Use this for an Nx1 strip.

If the strip geometry is 1xN, specifying three corners does not make sense. However, you should still specify three points. Corners 1, 2 and 3 should be the first three points starting from one endpoint of the line, such that corner 1 is at the endpoint, like such:

C1-C2-C3-...

If the 1xN strip has fewer than three electrodes, manual labeling should be used instead. You should also set the naming convention to line.

Following automatic labeling, it is recommended to inspect the automatically assigned geometry and make sure it is correct. Several common issues include rotation and reflection. To resolve these, select the rotation and reflection actions from the Labeling menu of the Electrode window.

It's important to take note of the labeling specific parameters at the bottom of the electrode localization window. For an explanation of what these parameters do, see below the section on [https://github.com/aestrivex/ielu/wiki#parameters-and-settings]. However, it's important because you might not want them by default to have the same values as the sorting parameters. Because we are more certain about the electrodes to be included in the algorithm, we probably want to reduce the value of delta, to avoid strange looking configurations. But if your grid has odd or aberrant curvature, like the one shown below, it might be better to increase delta.

Adding new points

See https://github.com/aestrivex/ielu/wiki#creating-new-electrodes above

ROI identification

You can identify what ROIs are contacted by specific electrodes. This procedure is unlikely to be as reliable as using some other algorithms, because the registration procedure has a potentially very high error.

Note - The ROIs are by default based on the Freesurfer segmentation for subcortical ROIs, and on the Freesurfer 'aparc' parcellation for cortical ROIs. However, you can instead use any other surface parcellation that you choose, as long as the annotations exist in SUBJECTS_DIR/SUBJECT/label/. To adjust the parcellation and the error radius, you can edit these fields in the menu resulting from ''Settings -> Visualizations/Output''.

The operation ''Tools -> Estimate all ROI contacts'' will do this for every electrode assigned to any grid.

Because of the registration error, this procedure doesn't really make a whole lot of sense, but some of my users really wanted it.

Registration error

The registration procedures we use do not always work perfectly. Here is one subject in which it failed:

The CT is shown in purple, while the MRI is shown in grayscale. If you look carefully, you can see that they do not align. The most prominent purple lines show the position of the skull, which is not aligned with the position of the skull in the MRI.

This view is shown with freeview, bundled with freesurfer. It shows the CT image registered to the MR space, and is overlaid on the MRI.

Typically the main problem with the pipeline is the registration between the CT and MR images. This cannot be fixed in all cases, and there is a little bit of an art to swapping the parameters, but we have some tools and procedures to get around some of the common problems we observed.

Fixing skewed registrations

One of the main problems with registration (especially on MGH data) is CT images which are skewed. That is, the dicom headers in the CT might indicate 1mm slice thickness, but the CT in actual fact appears to have a lower slice thickness. An example is shown below:

In the image above, notice that the shape of the brain (after the visualization program freeview accounts for the voxel sizes indicated in the image header) is quite skewed. The image below shows a CT with a correct shape.

A word of warning - Unfortunately, sometimes these routines are overzealous, and can take up to 15 minutes to run. Thus, it is important for the user to be able to make manual overrides.

Co-registration (CT to MR)

There are three possible registration settings in the options menu. By default, experimental shape correction is selected.

• Experimental shape correction
  • Assumes that the image is skewed along the Z-axis. X/Y information should be correct. It tries to use an algorithm to scale the image back to the correct shape in the Z direction.
  • If the user specifies "override_zoom_factor" as anything other than 0, the program takes this as an override. This causes the program to resample the image to the correct shape based on the zoom factor provided.
  • The zoom factor (zf) is an estimate of how skewed the image is. The image will be resampled by a factor of 1 / zf.
  • For the first image shown above, an approximately correct value of zf would be [1.3, 1, 1] (that is; the image is approximately 1.3 times larger than it should be in the Z direction, and no change is needed in the X and Y axes).
  • Note - the program automatically tries to determine what is the correct value of zf; the user is able to override the program's choice here in case the registration does not make sense.
• uncorrected MI registration
  • Does mutual information based registration without the resampling described above (and in the freesurfer documentation)
  • This is equivalent to manually setting a zoom factor of [1, 1, 1] (but saves the transformation separately from manually setting a zoom factor of 1. This is important because sometimes the algorithm overzealously tries to do zoom factor corrections when the zoom factor is 1, which takes time.)
• No registration
  • Assumes the CT has already been manually registered to the freesurfer space. (This may be buggy).

I recommend for the user to look at the CT images in freeview or fslview before running them and guesstimate what is the right zoom factor. A good guesstimate is usually good enough for the program to get it right.

Checking the automated registrations

I highly recommend examining the automated registrations by overlaying the automated registration and the T1 image on top of one another in freeview. This is a critical debugging tool, but also a good sanity check just for the user.

To do this you need to know where the registrations are. The registrations output by ielu are in the following locations:

• experimental shape correction : $SUBJECTS_DIR/$SUBJECT/mri/ct_reg/ct_final_resamp_reg.nii.gz
• uncorrected MI registration : $SUBJECTS_DIR/$SUBJECT/mri/ct_nas.nii.gz

(ct_nas is supposed to stand for Native Anatomical Space since the registration is done to rawavg and not orig).

You can check these registrations like so: freeview orig.mgz ct_nas.nii.gz

This will load both images but show one of them on top of the other. You should adjust the opacity of the one on top to make sure they align closely.

Snapping the ROIs to the surface

The algorithms to coregister the CT and MR images are imperfect, even when they work. In order to further correct for registration errors we utilize a snapping procedure described in Dykstra et al. (2012). This procedure is only effective on subdural electrodes placed on the surface, which can be localized as surface coordinates.

The snapping algorithm constructs imaginary springs between each pair of adjacent electrodes on the grid. These springs each have a preferred distance at which they like to exist; if the spring is stretched out too thinly then it exhibits a lot of potential energy, and if it is squished into a small distance it also contains a lot of potential energy. The snapping algorithm positions the electrodes along the cortical surface in such a way to search for the global minimum in a cost function that penalizes springs for having high potential energy, and penalizes electrodes from deviating too far from their starting location.

This procedure isn't perfect but works pretty well at reducing the registration error for some of the subdural grids.

In the electrode window, each grid has a type which is set to either "depth" or "subdural". This type is set to "depth" by default unless the grid went through the sorting and labeling algorithm in the pipeline; each grid to be snapped must be set to "subdural" in its electrode window.

To execute the snapping procedure, use the menu item Electrode -> Snap electrodes to surface in the main GUI. This will perform the snapping procedure on all of the electrodes listed as subdural grids.

Saving and exporting electrode locations

CSV Coordinates

Selecting Electrode -> Save CSV will save all electrodes that are assigned to a grid/strip. The output format is a CSV file with each row being an electrode, and columns being the electrode name, x, y, and z position.

Montages

Each grid is saved into its own separate montage file. The montage file is a plaintext ASCII file which contains the grid names and the electrode locations at that grid.

The montage file should be saved with the extension .sfp, which is not currently enforced by the software. However, if your montage file has an invalid file extension, MNE python will refuse to read it. This might be streamlined at some point later.

The following code describes how to convert a montage file into a .fif file containing the iEEG data

# read in the raw file. The raw files I have worked with were all in .edf format, but MNE python has
# readers for various formats
ra = mne.io.read_raw_edf( raw_edf_file )

# read the montage file
# the read_montage function has some very awkward syntax -- it can only read absolute paths or
# fully specified paths utilizing the path attribute. specifying
# a relative path will not work. This is because the path attribute has a default directory in the MNE
# python installation and is used to load different files containing device coordinates.
montage = mne.channels.read_montage( montage_file, path=directory_containing_montage_file)

# apply the electrode locations in the montage to the raw file
mne.channels.apply_montage( ra.info, montage )

# this function changes the channel type from EEG to sEEG in the mne data structure.
# this is not strictly necessary as most analysis software ignores this attribute.
# however it is important for loading data into the custom visualization tool inaivu because this tool allows
# in princple the combination of visualizing noninvasive scalp EEG and invasive electrodes simultaneously.
mne.channels.rename_channels( ra.info, dict((c, (c, 'seeg')) for c in montage.ch_names ))

# save the result to a fif file
ra.save( 'something.fif' )

Parameters and settings

Changing the default settings may be necessary. The settings window shows a lot of text trying to explain what the various parameters do, but here is some more detail and practical usage.

CT threshold

This controls the threshold used to extract electrodes from the CT image. I haven't figured out a good way to automatically determine this value across images -- some people suggested things like "3 standard deviations above the mean value of the CT image." Unfortunately, the CT images vary extensively in their contrast quality and consequently this type of formula didn't work at all. Consequently, I recommend to check the CT images manually, and pick a value that is lower than the maximum intensity of the electrodes, but higher than almost all of the skull (including lots of noise from the skull at the very least causes the visualization to become cluttered and will often cause the sorting procedure to fail).

Number of steps

This controls the number of iterations used in the snapping algorithm described in (Dykstra et al. 2012). In order to minimize the global cost function described in that paper, I use a simulated annealing algorithm with a specific (currently hard-coded) pretty rapid cooling schedule. This is a little different from the MATLAB implementation described in the paper uses a generic quadratic programming solver. Increasing the number of iterations should increase the robustness of the fit at the cost of the computation time.

Deformation constant

This constant controls the global cost function in the snapping algorithm described in (Dykstra et al. 2012). This alters the weight of the deformation term, which controls the absolute distance between the candidate snapped point and the starting point. Setting it >1 will reduce the possible range of vertices the point could take (by penalizing vertices further away more heavily). Setting this <1 will loosen this constraint.

Mask extracranial noise

Dilates the area of the brain (using a T1 brainmask, following registration of the electrodes to the T1 space) by a spherical kernel of radius 1 with the number of iterations specified in the settings panel. This dilated brainmask is used to mask noise very far from the brain.

Overwrite existing transformations

The transformation matrix between the CT and MR images for the given registration procedure are time consuming to generate and saved in a subdirectory of the freesurfer reconstruction. Checking this box will ensure to repeat the process and overwrite any old matrices.

Disable binary erosion procedure to reduce CT noise

The electrode locations are by default eroded with a spherical kernel of radius 1 (only once). This procedure is described in Taimouri et al. (2014). It is however completely undesirable with very low resolution CT images, and should be turned off.

Registration procedure

By default, CT and MR images are registered using an experimental algorithm to do the registration. This algorithm attempts to manually resample in the direction of the slice acquisition to correct for errors in head shape. First this does a number of intermediate registrations to figure out what is the correct factor to resample by, and then does the resampling, and registers the resampled image. The resampling procedure takes ten to fifteen minutes.

If set to regular MI registration, the resampling is skipped and the program will try to register the unadulterated CT image, which sometimes fails in our experience, particularly when the head shape is skewed and the CT image is far from isotropic.

When set to no registration, no registration is performed at all. This assumes a CT image which has already been manually registered to the subject's anatomical space.

Disable binary erosion

When unchecked, the program will automatically apply a binary erosion filter to remove noise from the CT image. This procedure only works with fairly high resolution CT images; If the image is low resolution (e.g., 256 voxels isotropic or lower), this will cause the entire image to be eroded. Check to disable when working with high resolution images such as those manually registered to the anatomical space.

Delta

This "delta" parameter is the fitting parameter that I found to be the most important. What it does is, when trying to automatically determine the local grid geometry, there are several constraints having to do with angles and distances to known candidate points. The distance constraint is that the distance d between candidate electrodes, must be as follows: f*(1-delta) < d < f*(1+delta), where f is the program's estimate of the fundamental distance between the electrodes on the individual grid being fit. Tweaking this parameter might provide better results, especially if you can see that the grid organization is either very ordered, or very erratic.

Epsilon

This parameter controls the threshold of how far away from 90 degrees is an acceptable starting point for the sorting and labeling algorithm. Increasing it allows the algorithm to run longer and possibly find a satisfactory solution.

Rho

This set of parameters control what is the acceptable deviation from the expected angles of a perfectly right angle NxM geometry grid. Changing these parameters doesn't tend to make a big difference.

Isotropization

By default the sorting and labeling algorithm operates and deals with electrodes in an isotropic coordinate space. Unfortunately, not all CT images are isotropic. Therefore, we translate them to an isotropic space:

• By header does the isotropization, by expanding the image to an isotropic coordinate space according to transformations derived from the coordinate space in the image header of the CT image
• By voxel assumes an isotropic image and makes it so that each dimension of the image has an identical number of voxels.
• Manual override allows the user to input a vector [xf, yf, zf] which creates an isotropic coordinate space by multiplying the axes dimensions by this vector. This vector will usually be the same as the registration override, if applicable.
• Isotropization off does not convert to an isotropic coordinate space.