{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# Corrupt known signal with point spread\n\nThe aim of this tutorial is to demonstrate how to put a known signal at a\ndesired location(s) in a :class:`mne.SourceEstimate` and then corrupt the\nsignal with point-spread by applying a forward and inverse solution.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Authors: The MNE-Python contributors.\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n\nimport mne\nfrom mne.datasets import sample\nfrom mne.minimum_norm import apply_inverse, read_inverse_operator\nfrom mne.simulation import simulate_evoked, simulate_stc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we set some parameters.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "seed = 42\n\n# parameters for inverse method\nmethod = \"sLORETA\"\nsnr = 3.0\nlambda2 = 1.0 / snr**2\n\n# signal simulation parameters\n# do not add extra noise to the known signals\nnave = np.inf\nT = 100\ntimes = np.linspace(0, 1, T)\ndt = times[1] - times[0]\n\n# Paths to MEG data\ndata_path = sample.data_path()\nsubjects_dir = data_path / \"subjects\"\nfname_fwd = data_path / \"MEG\" / \"sample\" / \"sample_audvis-meg-oct-6-fwd.fif\"\nfname_inv = data_path / \"MEG\" / \"sample\" / \"sample_audvis-meg-oct-6-meg-fixed-inv.fif\"\nfname_evoked = data_path / \"MEG\" / \"sample\" / \"sample_audvis-ave.fif\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load the MEG data\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fwd = mne.read_forward_solution(fname_fwd)\nfwd = mne.convert_forward_solution(fwd, force_fixed=True, surf_ori=True, use_cps=False)\nfwd[\"info\"][\"bads\"] = []\ninv_op = read_inverse_operator(fname_inv)\n\nraw = mne.io.read_raw_fif(data_path / \"MEG\" / \"sample\" / \"sample_audvis_raw.fif\")\nraw.info[\"bads\"] = []\nraw.set_eeg_reference(projection=True)\nevents = mne.find_events(raw)\nevent_id = {\"Auditory/Left\": 1, \"Auditory/Right\": 2}\nepochs = mne.Epochs(raw, events, event_id, baseline=(None, 0), preload=True)\nevoked = epochs.average()\n\nlabels = mne.read_labels_from_annot(\"sample\", subjects_dir=subjects_dir)\nlabel_names = [label.name for label in labels]\nn_labels = len(labels)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimate the background noise covariance from the baseline period\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "cov = mne.compute_covariance(epochs, tmin=None, tmax=0.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generate sinusoids in two spatially distant labels\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# The known signal is all zero-s off of the two labels of interest\nsignal = np.zeros((n_labels, T))\nidx = label_names.index(\"inferiorparietal-lh\")\nsignal[idx, :] = 1e-7 * np.sin(5 * 2 * np.pi * times)\nidx = label_names.index(\"rostralmiddlefrontal-rh\")\nsignal[idx, :] = 1e-7 * np.sin(7 * 2 * np.pi * times)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Find the center vertices in source space of each label\n\nWe want the known signal in each label to only be active at the center. We\ncreate a mask for each label that is 1 at the center vertex and 0 at all\nother vertices in the label. This mask is then used when simulating\nsource-space data.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "hemi_to_ind = {\"lh\": 0, \"rh\": 1}\nfor i, label in enumerate(labels):\n # The `center_of_mass` function needs labels to have values.\n labels[i].values.fill(1.0)\n\n # Restrict the eligible vertices to be those on the surface under\n # consideration and within the label.\n surf_vertices = fwd[\"src\"][hemi_to_ind[label.hemi]][\"vertno\"]\n restrict_verts = np.intersect1d(surf_vertices, label.vertices)\n com = labels[i].center_of_mass(\n subjects_dir=subjects_dir, restrict_vertices=restrict_verts, surf=\"white\"\n )\n\n # Convert the center of vertex index from surface vertex list to Label's\n # vertex list.\n cent_idx = np.where(label.vertices == com)[0][0]\n\n # Create a mask with 1 at center vertex and zeros elsewhere.\n labels[i].values.fill(0.0)\n labels[i].values[cent_idx] = 1.0\n\n # Print some useful information about this vertex and label\n if \"transversetemporal\" in label.name:\n dist, _ = label.distances_to_outside(subjects_dir=subjects_dir)\n dist = dist[cent_idx]\n area = label.compute_area(subjects_dir=subjects_dir)\n # convert to equivalent circular radius\n r = np.sqrt(area / np.pi)\n print(\n f\"{label.name} COM vertex is {dist * 1e3:0.1f} mm from edge \"\n f\"(label area equivalent to a circle with r={r * 1e3:0.1f} mm)\"\n )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create source-space data with known signals\n\nPut known signals onto surface vertices using the array of signals and\nthe label masks (stored in labels[i].values).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "stc_gen = simulate_stc(fwd[\"src\"], labels, signal, times[0], dt, value_fun=lambda x: x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot original signals\n\nNote that the original signals are highly concentrated (point) sources.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "kwargs = dict(\n subjects_dir=subjects_dir,\n hemi=\"split\",\n smoothing_steps=4,\n time_unit=\"s\",\n initial_time=0.05,\n size=1200,\n views=[\"lat\", \"med\"],\n)\nclim = dict(kind=\"value\", pos_lims=[1e-9, 1e-8, 1e-7])\nbrain_gen = stc_gen.plot(clim=clim, **kwargs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulate sensor-space signals\n\nUse the forward solution and add Gaussian noise to simulate sensor-space\n(evoked) data from the known source-space signals. The amount of noise is\ncontrolled by ``nave`` (higher values imply less noise).\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "evoked_gen = simulate_evoked(fwd, stc_gen, evoked.info, cov, nave, random_state=seed)\n\n# Map the simulated sensor-space data to source-space using the inverse\n# operator.\nstc_inv = apply_inverse(evoked_gen, inv_op, lambda2, method=method)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot the point-spread of corrupted signal\n\nNotice that after applying the forward- and inverse-operators to the known\npoint sources that the point sources have spread across the source-space.\nThis spread is due to the minimum norm solution so that the signal leaks to\nnearby vertices with similar orientations so that signal ends up crossing the\nsulci and gyri.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "brain_inv = stc_inv.plot(**kwargs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercises\n - Change the ``method`` parameter to either ``'dSPM'`` or ``'MNE'`` to\n explore the effect of the inverse method.\n - Try setting ``evoked_snr`` to a small, finite value, e.g. 3., to see the\n effect of noise.\n\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.2" } }, "nbformat": 4, "nbformat_minor": 0 }