Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add TFR vector visualization #6

Open
wants to merge 7 commits into
base: main
Choose a base branch
from
Open

Add TFR vector visualization #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
c9992fd
Add TFR vector visualization
alexrockhill Apr 3, 2023
f2d105e
fix style, try to fix 3D rendering
alexrockhill Apr 3, 2023
1b70bc8
MAINT: Make implicit req explicit
larsoner Apr 4, 2023
9fa75a7
FIX: Black
larsoner Apr 4, 2023
afe2c34
rerun cis
alexrockhill Apr 5, 2023
213016a
FIX: CIs
larsoner Apr 5, 2023
e8de33f
FIX: CIs
larsoner Apr 5, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .circleci/config.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ jobs:
echo "set -e" >> $BASH_ENV
echo "export OPENBLAS_NUM_THREADS=4" >> $BASH_ENV
echo "export XDG_RUNTIME_DIR=/tmp/runtime-circleci" >> $BASH_ENV
echo "export MNE_3D_OPTION_MULTI_SAMPLES=1" >> $BASH_ENV
echo "export PATH=~/.local/bin/:$PATH" >> $BASH_ENV
echo "export DISPLAY=:99" >> $BASH_ENV
echo "source ~/python_env/bin/activate" >> $BASH_ENV
Expand Down
8 changes: 7 additions & 1 deletion .github/workflows/tests.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,13 @@ jobs:
if: "matrix.mne-version == 'mne-main'"
run: git clone --single-branch --branch main https://github.com/mne-tools/mne-python.git
- run: pip install -ve ./mne-python
- run: pip install -v ${{ matrix.qt }}
- run: |
QT_LIB=${{ matrix.qt }}
if [[ "${{ matrix.qt }}" == "PySide6" ]]; then
QT_LIB="PySide6!=6.5"
fi
pip install -v "$QT_LIB"
shell: bash -e {0}
- run: pip install -ve .[tests]
- run: mne sys_info
- run: |
Expand Down
31 changes: 26 additions & 5 deletions mne_gui_addons/_vol_stc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,14 @@ def _int_complex_conj(data):
return (conj // (conj.max() // RANGE_VALUE + 1)).astype(BASE_INT_DTYPE)


def _complex_to_vector(data):
if data.dtype == COMPLEX_DTYPE:
data = data["re"] + 1j * data["im"]
assert np.iscomplexobj(data)
data = np.abs(data) * np.cos(np.angle(data))
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 equivalent to data.real

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> data = rng.normal(size=(1000,)) + 1j * rng.normal(size=(1000,))
>>> np.allclose(data.real, np.abs(data) * np.cos(np.angle(data)))
True

Copy link
Collaborator Author

@alexrockhill alexrockhill Apr 6, 2023 •
edited
Loading

Choose a reason for hiding this comment

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

Interesting, so that was just plotting the real component, that's even simpler to understand potentially. What do you think about interpreting this? The real component is just the amplitude of the cosine for that frequency compared to what I was originally thinking was tapering the magnitude of the vector by it's alignment with 0 or 180 degree phase. Since those are equivalent, I'm thinking there might be a simpler interpretation, I guess that the cosine is already zero at 90 and 270 degrees so the real part already includes that tapering.

Copy link
Member

Choose a reason for hiding this comment

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

I have a hard time interpreting the real or imaginary parts of these filters. In spectral analysis it means some cosine or sine part of a complex exponential (it's just a convenient way of bookkeeping). And in cross-spectrum / coherence analysis it means some phase shift between signals of interest. So plotting just the real (or imag) part of the signal only gives you part of the information. I don't know the right thing to do here but plotting just the real part (or just the imaginary part) seems to miss a lot

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Hmm those are good examples of where the real and imaginary parts of the filters are used directly but I think via the Hilbert transform, the imaginary part of the signal is just the real signal shifted in phase by a quarter cycle so I think that way the information is actually redundant and you might not be missing anything after all. In this gist (https://gist.github.com/alexrockhill/58c4316415c54fcc24f85a617862e2bf) I simulate a pure sine wave and then plot it via this method and similarly in the code below, you what I would expect to recover; a vector oscillating in the direction that was simulated. If there is a reason why visualizing it this way and this interpretation is wrong, I'm happy to abandon the idea or table until we do it some other way that is better, but this is what I was thinking would be helpful to visualize--when you get a pure oscillation, the vector changes direction forward and backward on the brain in the direction of that oscillation, I think the math on how to get there changes the interpretation but I don't quite understand why this isn't helpful or what is wrong about this visualization.

import numpy as np
import mne
import matplotlib.pyplot as plt
info = mne.create_info(['1'], 1000, 'eeg')
epochs = mne.EpochsArray(np.sin(np.linspace(0, 10, 10001) * 2 * np.pi * 20)[None, None] * 1e-6, info)
epochs_tfr = mne.time_frequency.tfr_morlet(epochs, [20], [10], return_itc=False, average=False, output='complex')
plt.plot(epochs_tfr.data.real.flatten())
plt.plot(epochs_tfr.data.imag.flatten())

Copy link
Member

Choose a reason for hiding this comment

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

the imaginary part of the signal is just the real signal shifted in phase by a quarter cycle

Agreed. Same can be said of the real and imaginary parts of the Fourier basis (cos +j sin).

I think that way the information is actually redundant

No, I don't think this is how the Hilbert transform is meant to be interpreted/thought of in practice in signal processing any more than you can think of the real and imaginary coefficients of a Fourier decomposition as being redundant (the bases are actually orthogonal). And it becomes more complicated once you're talking about it in the context of a cross special density estimate, which underlies DICS, which is where I assume this would actually be useful. I think a deeper understanding of the underlying signal processing principles at play here (ideally derived from math) should guide any public implementation here rather than simulations.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

The imaginary part isn't totally redundant, but it is similar to the phase-shifted real component in actual data. The cosine of the phase scaling the magnitude is identical to the real value in the experimental data also.

import numpy as np
import mne
import matplotlib.pyplot as plt
evoked = mne.read_evokeds(mne.datasets.sample.data_path() /  'MEG' / 'sample' / 'sample_audvis-ave.fif')[0]
evoked_tfr = mne.time_frequency.tfr_morlet(evoked, [20], [2], return_itc=False, average=False, output='complex')
n = evoked.info['sfreq'] / 20 / 4
plt.plot(evoked_tfr.data[0, 0, 0].real)
plt.plot(evoked_tfr.data[0, 0, 0, int(round(n)):].imag)
plt.show()

np.allclose(evoked_tfr.data.real, np.abs(evoked_tfr.data) * np.cos(np.angle(evoked_tfr.data)))

As I understand it, the cross-spectral density is just for whitening and the complex data itself is what gets input to the minimum norm or beamformer solutions.

I don't want to post hoc justify that the real component is what is the right thing to visualize when the phase-scaled magnitude was what I was going for and was theoretically motivated but it would also be nice to visualize the direction. I know a first-component-of-the-SVD method has also been proposed which would be interesting to see implemented but it would also be nice to plot something that isn't just the direction of maximum variance if possible. Since this gets a bit complicated to understand and interpret, it would be nice if there was something with minimal processing of the data. The real component is as minimal as it gets but obviously might be too simple and not the right solution.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Thanks for a good discussion, how about we table it for now and come back to it later?

Copy link
Member

Choose a reason for hiding this comment

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

Yeah we can discuss it later

return data


class VolSourceEstimateViewer(SliceBrowser):
"""View a source estimate time-course time-frequency visualization."""

Expand Down Expand Up @@ -326,9 +334,18 @@ def __init__(
)

self._data_max = abs(stc_data).max()
if self._data.shape[2] > 1 and not self._is_complex:
if self._data.shape[2] > 1:
# also compute vectors for chosen time
self._stc_vectors = self._pick_stc_tfr(stc_data).astype(float)
if self._f_idx is None:
self._stc_vectors = self._pick_stc_tfr(data)
self._tfr_vector_max = None
else:
self._stc_vectors = _complex_to_vector(data)
self._tfr_vector_max = self._stc_vectors.max()
self._stc_vectors = self._pick_stc_tfr(self._stc_vectors)
self._stc_vectors = self._pick_epoch(self._stc_vectors).astype(float)
if self._f_idx is not None:
self._stc_vectors /= self._tfr_vector_max
self._stc_vectors /= self._data_max
self._stc_vectors_masked = self._stc_vectors.copy()

Expand Down Expand Up @@ -428,10 +445,14 @@ def _update_stc_pick(self):
self._stc_range = max([stc_max, -self._stc_min]) - self._stc_min

def _update_vectors(self):
if self._data.shape[2] > 1 and not self._is_complex:
if self._data.shape[2] > 1:
# pick vector as well
self._stc_vectors = self._pick_stc_tfr(self._data)
if self._f_idx is not None:
self._stc_vectors = _complex_to_vector(self._stc_vectors)
self._stc_vectors = self._pick_epoch(self._stc_vectors).astype(float)
if self._f_idx is not None:
self._stc_vectors /= self._tfr_vector_max
self._stc_vectors /= self._data_max
self._update_vector_threshold()
self._plot_vectors()
Expand Down Expand Up @@ -1236,7 +1257,7 @@ def _update_cmap(
if not update_3d:
return

if self._data.shape[2] > 1 and not self._is_complex:
if self._data.shape[2] > 1:
# update vector mask
self._update_vector_threshold()
self._plot_vectors(draw=False)
Expand Down Expand Up @@ -1343,7 +1364,7 @@ def set_3d_view(

def _plot_vectors(self, draw=True):
"""Update the vector plots."""
if self._data.shape[2] > 1 and not self._is_complex:
if self._data.shape[2] > 1:
self._vector_data.point_data["vec"] = (
VECTOR_SCALAR * self._stc_vectors_masked
)
Expand Down
1 change: 1 addition & 0 deletions pyproject.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ dependencies = [
"qtpy",
"pyvista",
"pyvistaqt",
"matplotlib",
"mne",
"nibabel",
"dipy>=1.4",
Expand Down