Repository: mne-tools/mne-nirs (language: python)

Issue:
Add ability to interface stc_to_sensor with GLM results



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/CircleCI_PR.yml b/.github/workflows/CircleCI_PR.yml
index 9789249c4..222c3724d 100644
--- a/.github/workflows/CircleCI_PR.yml
+++ b/.github/workflows/CircleCI_PR.yml
@@ -8,5 +8,5 @@ jobs:
         uses: larsoner/circleci-artifacts-redirector-action@master
         with:
           repo-token: ${{ secrets.GITHUB_TOKEN }}
-          artifact-path: 0/dev/master/index.html
+          artifact-path: 0/dev/pull/${{ github.event.pull_request.number }}/index.html
           circleci-jobs: build_docs
diff --git a/doc/api.rst b/doc/api.rst
index 354a5b63d..50106400f 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -91,6 +91,7 @@ Visualisation
    plot_glm_contrast_topo
    plot_glm_group_topo
    plot_timechannel_quality_metric
+   plot_glm_surface_projection
 
 
 Simulation
diff --git a/doc/conf.py b/doc/conf.py
index 0aea7d846..fb083bc00 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -225,6 +225,7 @@
 sphinx_gallery_conf = {
     'doc_module': 'mne_nirs',
     'backreferences_dir': os.path.join('generated'),
+    'image_scrapers': ('pyvista', 'matplotlib'),
     'reference_url': {
         'mne_nirs': None},
     'download_all_examples': False,
diff --git a/environment.yml b/environment.yml
index 7c1282601..436b8a322 100644
--- a/environment.yml
+++ b/environment.yml
@@ -33,10 +33,10 @@ dependencies:
 - nibabel
 - python-picard
 - pyqt!=5.15.3
-- mne>=0.21.0
 - mffpy>=0.5.7
 - ipywidgets
 - lxml
 - pip:
   - ipyvtklink
   - https://github.com/nilearn/nilearn/archive/main.zip
+  - https://github.com/mne-tools/mne-python/archive/main.zip
diff --git a/examples/general/plot_12_group_glm.py b/examples/general/plot_12_group_glm.py
index 25f6b1aba..91412c634 100644
--- a/examples/general/plot_12_group_glm.py
+++ b/examples/general/plot_12_group_glm.py
@@ -89,6 +89,7 @@
 from mne_nirs.utils._io import glm_to_tidy
 from mne_nirs.visualisation import plot_glm_group_topo
 from mne_nirs.datasets import fnirs_motor_group
+from mne_nirs.visualisation import plot_glm_surface_projection
 
 # Import MNE-BIDS processing
 from mne_bids import BIDSPath, read_raw_bids
@@ -393,6 +394,56 @@ def individual_analysis(bids_path, ID):
                     con_model_df, colorbar=True, threshold=True)
 
 
+###############################################################################
+# Cortical Surface Projections
+# ----------------------------
+#
+# .. sidebar:: 3D plotting
+#
+#    Due to limitations with the cloud computing used to generate these
+#    documents, the 3D plot will not appear in the online documentation.
+#    However, when you run this code locally you will obtain a 3D head figure
+#    simillar to what is shown in the tutorial:
+#    `ECOG tutorial <https://mne.tools/dev/auto_tutorials/clinical/30_ecog.html#sphx-glr-auto-tutorials-clinical-30-ecog-py>`_.
+#
+# The topographic plots above can sometimes be difficult to interpret with
+# respect to the underlying cortical locations. It is also possible to present
+# the data by projecting the channel level GLM values to the nearest cortical
+# surface. This can make it easier to understand the spatial aspects of your
+# data. Note however, that this is not a complete forward model with photon
+# migration simulations.
+# In the figure below we project the group results from the two conditions
+# to the cortical surface, and also present the contrast results in the same
+# fashion.
+# As in the topo plots above you can see that the activity is predominately
+# contralateral to the side of finger tapping.
+
+
+# Generate brain figure from data
+clim = dict(kind='value', pos_lims=(0, 8, 11))
+brain = plot_glm_surface_projection(raw_haemo.copy().pick("hbo"),
+                                    con_model_df, clim=clim, view='dorsal',
+                                    colorbar=True, size=(800, 700))
+brain.add_text(0.05, 0.95, "Left-Right", 'title', font_size=16, color='k')
+
+# Run model code as above
+clim = dict(kind='value', pos_lims=(0, 11.5, 17))
+for idx, cond in enumerate(['Tapping/Left', 'Tapping/Right']):
+
+    # Run same model as explained in the sections above
+    ch_summary = df_cha.query("Condition in [@cond]")
+    ch_summary = ch_summary.query("Chroma in ['hbo']")
+    ch_model = smf.mixedlm("theta ~ -1 + ch_name", ch_summary,
+                           groups=ch_summary["ID"]).fit(method='nm')
+    model_df = statsmodels_to_results(ch_model)
+
+    # Generate brain figure from data
+    brain = plot_glm_surface_projection(raw_haemo.copy().pick("hbo"),
+                                        model_df, clim=clim, view='dorsal',
+                                        colorbar=True, size=(800, 700))
+    brain.add_text(0.05, 0.95, cond, 'title', font_size=16, color='k')
+
+
 ###############################################################################
 # Table of channel level results
 # ------------------------------
diff --git a/mne_nirs/visualisation/__init__.py b/mne_nirs/visualisation/__init__.py
index 936805e57..5629868c5 100644
--- a/mne_nirs/visualisation/__init__.py
+++ b/mne_nirs/visualisation/__init__.py
@@ -7,4 +7,5 @@
 from ._plot_nirs_source_detector import plot_nirs_source_detector
 from._plot_GLM_topo import (plot_glm_topo, plot_glm_contrast_topo,
                             plot_glm_group_topo)
-from ._plot_quality_metrics import plot_timechannel_quality_metric
\ No newline at end of file
+from ._plot_quality_metrics import plot_timechannel_quality_metric
+from ._plot_GLM_surface_projection import plot_glm_surface_projection
\ No newline at end of file
diff --git a/mne_nirs/visualisation/_plot_GLM_surface_projection.py b/mne_nirs/visualisation/_plot_GLM_surface_projection.py
new file mode 100644
index 000000000..788a4dba8
--- /dev/null
+++ b/mne_nirs/visualisation/_plot_GLM_surface_projection.py
@@ -0,0 +1,128 @@
+# Authors: Robert Luke <mail@robertluke.net>
+#
+# License: BSD (3-clause)
+
+import os
+import numpy as np
+
+from mne import stc_near_sensors, EvokedArray, read_source_spaces
+from mne.datasets import sample
+from mne.io.constants import FIFF
+from mne import verbose
+
+
+@verbose
+def plot_glm_surface_projection(raw, statsmodel_df, picks="hbo",
+                                value="Coef.",
+                                background='w', figure=None, clim='auto',
+                                mode='weighted', colormap='RdBu_r',
+                                surface='pial', hemi='both', size=800,
+                                view=None, colorbar=True, distance=0.03,
+                                subjects_dir=None, src=None, verbose=False):
+    """
+    Project GLM results on to the surface of the brain.
+
+    Note: This function provides a convenient wrapper around low level
+    MNE-Python functions. It is convenient if you wish to use a generic head
+    model. If you have acquired fMRI images you may wish to use the underlying
+    lower level functions.
+
+    Note: This function does not conduct a forward model analysis with photon
+    migration etc. It simply projects the values from each channel to the
+    local cortical surface. It is useful for visualisation, but users should
+    take this in to consideration when drawing conclusions from the
+    visualisation.
+
+    Parameters
+    ----------
+    raw : instance of Raw
+        Haemoglobin data.
+    statsmodel_df : dataframe
+        As produced by produced by `statsmodels_to_results`.
+    %(picks_base)s good sEEG, ECoG, and DBS channels.
+    value : str
+        Column from dataframe to plot.
+    background : matplotlib color
+        Color of the background of the display window.
+    figure : mayavi.core.api.Scene, matplotlib.figure.Figure, list, int, None
+        If None, a new figure will be created. If multiple views or a
+        split view is requested, this must be a list of the appropriate
+        length. If int is provided it will be used to identify the Mayavi
+        figure by it's id or create a new figure with the given id. If an
+        instance of matplotlib figure, mpl backend is used for plotting.
+    %(clim)s
+    mode : str
+        Can be "sum" to do a linear sum of weights, "weighted" to make this
+        a weighted sum, "nearest" to
+        use only the weight of the nearest sensor, or "single" to
+        do a distance-weight of the nearest sensor. Default is "sum".
+    colormap : str
+        Colormap to use.
+    surface : str
+        The type of surface (inflated, white etc.).
+    hemi : str
+        Hemisphere id (ie 'lh', 'rh', 'both', or 'split'). In the case
+        of 'both', both hemispheres are shown in the same window.
+        In the case of 'split' hemispheres are displayed side-by-side
+        in different viewing panes.
+    size : float or tuple of float
+        The size of the window, in pixels. can be one number to specify
+        a square window, or the (width, height) of a rectangular window.
+        Has no effect with mpl backend.
+    view : str
+        View to set brain to.
+    colorbar : bool
+        If True, display colorbar on scene.
+    distance : float
+        Distance (m) defining the activation "ball" of the sensor.
+    %(subjects_dir)s
+    src : instance of SourceSpaces
+        The source space.
+    %(verbose)s
+
+    Returns
+    -------
+    figure : instance of mne.viz.Brain | matplotlib.figure.Figure
+        An instance of :class:`mne.viz.Brain` or matplotlib figure.
+    """
+    if not (raw.ch_names == list(statsmodel_df['ch_name'].values)):
+        raise RuntimeError('MNE data structure does not match dataframe '
+                           f'results. Raw = {raw.ch_names}. '
+                           f'GLM = {list(statsmodel_df["ch_name"].values)}')
+
+    ea = EvokedArray(np.tile(statsmodel_df[value].values.T, (1, 1)).T,
+                     raw.copy().info)
+
+    # TODO: mimic behaviour of other MNE-NIRS glm plotting options
+    if picks is not None:
+        ea = ea.pick(picks=picks)
+
+    if subjects_dir is None:
+        subjects_dir = sample.data_path() + '/subjects'
+    if src is None:
+        fname_src_fs = os.path.join(subjects_dir, 'fsaverage', 'bem',
+                                    'fsaverage-ico-5-src.fif')
+        src = read_source_spaces(fname_src_fs)
+
+    picks = np.arange(len(ea.info['ch_names']))
+
+    # Set coord frame
+    for idx in range(len(ea.ch_names)):
+        ea.info['chs'][idx]['coord_frame'] = FIFF.FIFFV_COORD_HEAD
+
+    # Generate source estimate
+    kwargs = dict(
+        evoked=ea, subject='fsaverage', trans='fsaverage',
+        distance=distance, mode=mode,
+        subjects_dir=subjects_dir, src=src, project=True)
+    stc = stc_near_sensors(picks=picks, **kwargs, verbose=verbose)
+
+    # Produce brain plot
+    brain = stc.plot(src=src, subjects_dir=subjects_dir, hemi=hemi,
+                     surface=surface, initial_time=0, clim=clim, size=size,
+                     colormap=colormap, figure=figure, background=background,
+                     colorbar=colorbar, verbose=verbose)
+    if view is not None:
+        brain.show_view(view)
+
+    return brain
diff --git a/requirements.txt b/requirements.txt
index 3c19809ad..33ea6b69d 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -35,5 +35,5 @@ ipywidgets
 ipyvtklink
 
 # mne-nirs specific requirements
-mne>=0.21.0
+https://github.com/mne-tools/mne-python/archive/main.zip
 lxml
diff --git a/tools/github_actions_dependencies.sh b/tools/github_actions_dependencies.sh
index c277ecd64..7b9414d44 100755
--- a/tools/github_actions_dependencies.sh
+++ b/tools/github_actions_dependencies.sh
@@ -22,6 +22,7 @@ else # pip 3.9 (missing statsmodels and dipy)
 	python -c "import pyvistaqt"
 	pip install --progress-bar off --upgrade -r requirements.txt
 fi
+pip install --progress-bar off --upgrade https://github.com/mne-tools/mne-python/archive/main.zip
 pip install --progress-bar off --upgrade -r requirements_testing.txt
 if [ "${DEPS}" != "minimal" ]; then
 	pip install nitime


Test patch:
diff --git a/mne_nirs/visualisation/tests/test_visualisation.py b/mne_nirs/visualisation/tests/test_visualisation.py
index c4c55891c..ef03c12e8 100644
--- a/mne_nirs/visualisation/tests/test_visualisation.py
+++ b/mne_nirs/visualisation/tests/test_visualisation.py
@@ -11,7 +11,8 @@
     _load_dataset
 from mne_nirs.experimental_design import make_first_level_design_matrix
 from mne_nirs.statistics import run_GLM
-from mne_nirs.visualisation import plot_glm_topo
+from mne_nirs.visualisation import plot_glm_topo, plot_glm_surface_projection
+from mne_nirs.utils import glm_to_tidy
 
 
 @requires_pysurfer
@@ -161,3 +162,28 @@ def test_fig_from_axes():
     from mne_nirs.visualisation._plot_GLM_topo import _get_fig_from_axes
     with pytest.raises(RuntimeError, match="Unable to extract figure"):
         _get_fig_from_axes([1, 2, 3])
+
+
+@requires_pysurfer
+@traits_test
+@pytest.mark.filterwarnings('ignore:.*nilearn.glm module is experimental.*:')
+def test_run_plot_GLM_projection():
+    mne.viz.set_3d_backend('pyvista')
+    raw_intensity = _load_dataset()
+    raw_intensity.crop(450, 600)  # Keep the test fast
+
+    design_matrix = make_first_level_design_matrix(raw_intensity,
+                                                   drift_order=1,
+                                                   drift_model='polynomial')
+    raw_od = mne.preprocessing.nirs.optical_density(raw_intensity)
+    raw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od)
+    glm_estimates = run_GLM(raw_haemo, design_matrix)
+    df = glm_to_tidy(raw_haemo, glm_estimates, design_matrix)
+    df = df.query("Chroma in 'hbo'")
+    df = df.query("Condition in 'A'")
+
+    brain = plot_glm_surface_projection(raw_haemo.copy().pick("hbo"),
+                                        df, clim='auto', view='dorsal',
+                                        colorbar=True, size=(800, 700),
+                                        value="theta")
+    assert type(brain) == mne.viz._brain.Brain