Skip to content

Conversation

@mscheltienne
Copy link
Member

@mscheltienne mscheltienne commented Aug 25, 2022

Closes #11079
Closes #11088

On main, interpolate_bridged_electrodes creates one virtual channel for each pair of bridged electrodes. On this PR, it creates one virtual channel for each group of bridged electrodes. A warning is emitted if the group includes more than 4 electrodes (not included).

On main, the virtual channel location is determined by averaging the 3 spherical coordinates of each pair of electrodes. On this PR, the location is determined by averaging the 3 cartesian coordinates of each electrode in the group. The resulting position is converted to spherical coordinates and the radius is replaced with the average radius of the electrodes in the group before being converted back to cartesian coordinates.

I ran this example locally, just in case: https://mne.tools/stable/auto_examples/preprocessing/eeg_bridging.html#ex-eeg-bridging

@mscheltienne mscheltienne marked this pull request as draft August 25, 2022 19:16
Comment on lines 118 to 134
# look for group of bad channels
nodes = sorted(set(chain(*bridged_idx)))
G_dense = np.zeros((len(nodes), len(nodes)))
# fill the edges with a weight of 1
for bridge in bridged_idx:
idx0 = nodes.index(bridge[0])
idx1 = nodes.index(bridge[1])
G_dense[idx0, idx1] = 1
G_dense[idx1, idx0] = 1
# look for connected components
_, labels = connected_components(csr_matrix(G_dense), directed=False)
groups_idx = [
[nodes[j] for j in np.where(labels == k)[0]] for k in set(labels)
]
groups_names = [
[inst.info.ch_names[k] for k in group_idx] for group_idx in groups_idx
]
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This snippet looks for group of electrodes. It's equivalent to the networkx code:

G = nx.Graph()
for bridge in bridged_idx:
    G.add_edge(*bridge)
groups_idx = nx.connected_components(G)

I had no idea scipy had those functionalities, and I never used them. Maybe there is a nicer way to do this?
Also, I used connected_components(csr_matrix(G_dense), directed=False) as per the example of connected_components, but it also works without csr_matrix. Not sure if it's needed..

Comment on lines 191 to 226
def _find_centroid(ch_pos, group_names):
"""Compute the centroid position between N electrodes.
The centroid should be determined in spherical "head" coordinates which is
more accurante than cutting through the scalp by averaging in cartesian
coordinates.
A simple way is to average the location in cartesian coordinate, convert
to spehrical coordinate and replace the radius with the average radius of
the N points in spherical coordinates.
Parameters
----------
ch_pos : OrderedDict
The position of all channels in cartesian coordinates.
group_names :
The name of the N electrodes used to determine the centroid.
Returns
-------
pos_centroid : array of shape (3,)
The position of the centroid in cartesian coordinates.
"""
cartesian_positions = np.zeros((len(group_names), 3))
sphere_positions = np.zeros((len(group_names), 3))
for i, ch_name in enumerate(group_names):
cartesian_positions[i, :] = ch_pos[ch_name]
sphere_positions[i, :] = _cart_to_sph(ch_pos[ch_name])
cartesian_pos_centroid = np.average(cartesian_positions, axis=0)
sphere_pos_centroid = _cart_to_sph(cartesian_pos_centroid)
# average the radius and overwrite it
avg_radius = np.average(sphere_positions, axis=0)[0]
sphere_pos_centroid[0, 0] = avg_radius
# convert back to cartesian
pos_centroid = _sph_to_cart(sphere_pos_centroid)[0, :]
return pos_centroid
Copy link
Member Author

@mscheltienne mscheltienne Aug 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I moved the computation of an average location in a separate function because I wanted to add a test specifically for this. Maybe this function and the associated test could be moved to a different module? Not sure where would be a good fit..

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is an okay place for it, we can move if it's useful elsewhere

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't mind, but we will definitely forget it 😄

@mscheltienne mscheltienne changed the title [WIP] Improve interpolation of bridged electrodes [MRG] Improve interpolation of bridged electrodes Aug 26, 2022
@mscheltienne mscheltienne marked this pull request as ready for review August 26, 2022 11:09
Comment on lines 120 to 128
G_dense = np.zeros((len(nodes), len(nodes)))
# fill the edges with a weight of 1
for bridge in bridged_idx:
idx0 = nodes.index(bridge[0])
idx1 = nodes.index(bridge[1])
G_dense[idx0, idx1] = 1
G_dense[idx1, idx0] = 1
# look for connected components
_, labels = connected_components(csr_matrix(G_dense), directed=False)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would do this with https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.coo_matrix.html

Or if you know the matrices will always be small (they probably will) then what you have is fine. I wouldn't bother casting explicitly to csr_matrix, though

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It will always be small, a worst (impossible) case scenario would be a 256 electrodes cap with all electrodes bridged. So 256x256.

So coo_matrix to avoid the for loop and remove the csr_matrix. I'll check it later today.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Up to you if you want to coo_matrix or not. It looks like you might have to loop to get idx0 and idx1. In principle you could maybe use np.searchsorted to avoid even those loops but it starts to be less and less readable...

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll stick with the current version, I just changed the .index with a better np.searchsorted. I did not find an easy and readable way to convert the information to the format for coo_matrix. I was starting with:

nodes = sorted(set(chain(*bridged_idx)))
rows, cols = zip(*bridged_idx)
# need to convert the idx as with `nodes.index(idx)`
coo_matrix(np.ones(len(rows)), rows, cols, shape=(len(nodes), len(nodes)))
# ...

It doesn't look more readable, and the conversion of the idx requires the same loop.

@mscheltienne
Copy link
Member Author

Any way to kill linux / pip-pre / linux pip 3.10? It's been running for 5 hours..
How about adding timeout-minutes: 120 here

job:
name: 'linux pip 3.10'
runs-on: ubuntu-20.04

Copy link
Member

@larsoner larsoner left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just minor comments from me, otherwise LGTM!

@mscheltienne
Copy link
Member Author

@larsoner I addressed your last comments, good to merge on my end.
Do you want to wait for another review, maybe @alexrockhill as he authored the original function?

@larsoner
Copy link
Member

@mscheltienne I disabled the GitHub "notify on commits pushed" so I didn't see that you had addressed my comments (too many emails) so thanks for the ping! I think it's fine to merge as-is, but @alexrockhill let me know if you want to look

@larsoner larsoner enabled auto-merge (squash) August 29, 2022 21:33
@larsoner larsoner merged commit 4787bac into mne-tools:main Aug 29, 2022
@mscheltienne mscheltienne deleted the interpolate-bridged branch August 30, 2022 09:53
larsoner added a commit to larsoner/mne-python that referenced this pull request Aug 30, 2022
* upstream/main: (25 commits)
  DOC: Exclude some implicit refs in doc build (mne-tools#10433)
  MAINT: Better issue forms (mne-tools#11101)
  [FIX] Typo in example (mne-tools#11118)
  [MRG] Improve interpolation of bridged electrodes (mne-tools#11094)
  BUG: Spectrum deprecation cleanup [circle deploy] (mne-tools#11115)
  Add API entry list and map (mne-tools#10999)
  Add legacy decorator (mne-tools#11097)
  [ENH, MRG] Add time-frequency epoch source estimation (mne-tools#11095)
  Revert "Add error message when conversion of EEG locs to [circle deploy] (mne-tools#11104)
  MRG: Fixes for mne-tools#11090 (mne-tools#11108)
  add test for edf units param (mne-tools#11105)
  BUG: Improve logic for bti (mne-tools#11102)
  add spectrum class (mne-tools#10184)
  ENH : add units parameter to read_raw_edf in case units is missing from the file (mne-tools#11099)
  ENH: Add temperature and galvanic (mne-tools#11090)
  Add error message when conversion of EEG locs to head space fails (mne-tools#11080)
  DOC: removed unnecessary line in PSF example (mne-tools#11085)
  FIX: Update helmet during ICP (mne-tools#11084)
  Fix various typos (mne-tools#11086)
  DOC: Rel
  ...
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

3 participants