{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# Background on projectors and projections\n\nThis tutorial provides background information on projectors and Signal Space\nProjection (SSP), and covers loading and saving projectors, adding and removing\nprojectors from Raw objects, the difference between \"applied\" and \"unapplied\"\nprojectors, and at what stages MNE-Python applies projectors automatically.\n\nWe'll start by importing the Python modules we need; we'll also define a short\nfunction to make it easier to make several plots that look similar:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import os\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D # noqa\nfrom scipy.linalg import svd\nimport mne\n\n\ndef setup_3d_axes():\n ax = plt.axes(projection='3d')\n ax.view_init(azim=-105, elev=20)\n ax.set_xlabel('x')\n ax.set_ylabel('y')\n ax.set_zlabel('z')\n ax.set_xlim(-1, 5)\n ax.set_ylim(-1, 5)\n ax.set_zlim(0, 5)\n return ax" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What is a projection?\n\nIn the most basic terms, a *projection* is an operation that converts one set\nof points into another set of points, where repeating the projection\noperation on the resulting points has no effect. To give a simple geometric\nexample, imagine the point $(3, 2, 5)$ in 3-dimensional space. A\nprojection of that point onto the $x, y$ plane looks a lot like a\nshadow cast by that point if the sun were directly above it:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ax = setup_3d_axes()\n\n# plot the vector (3, 2, 5)\norigin = np.zeros((3, 1))\npoint = np.array([[3, 2, 5]]).T\nvector = np.hstack([origin, point])\nax.plot(*vector, color='k')\nax.plot(*point, color='k', marker='o')\n\n# project the vector onto the x,y plane and plot it\nxy_projection_matrix = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 0]])\nprojected_point = xy_projection_matrix @ point\nprojected_vector = xy_projection_matrix @ vector\nax.plot(*projected_vector, color='C0')\nax.plot(*projected_point, color='C0', marker='o')\n\n# add dashed arrow showing projection\narrow_coords = np.concatenate([point, projected_point - point]).flatten()\nax.quiver3D(*arrow_coords, length=0.96, arrow_length_ratio=0.1, color='C1',\n linewidth=1, linestyle='dashed')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

#### Note

The @ symbol indicates matrix multiplication on NumPy arrays, and was\n introduced in Python 3.5 / NumPy 1.10. The notation plot(*point) uses\n Python argument expansion_ to \"unpack\" the elements of point into\n separate positional arguments to the function. In other words,\n plot(*point) expands to plot(3, 2, 5).

\n\nNotice that we used matrix multiplication to compute the projection of our\npoint $(3, 2, 5)$onto the $x, y$ plane:\n\n\\begin{align}\\left[\n \\begin{matrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & 0 & 0 \\end{matrix}\n \\right]\n \\left[ \\begin{matrix} 3 \\\\ 2 \\\\ 5 \\end{matrix} \\right] =\n \\left[ \\begin{matrix} 3 \\\\ 2 \\\\ 0 \\end{matrix} \\right]\\end{align}\n\n...and that applying the projection again to the result just gives back the\nresult again:\n\n\\begin{align}\\left[\n \\begin{matrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & 0 & 0 \\end{matrix}\n \\right]\n \\left[ \\begin{matrix} 3 \\\\ 2 \\\\ 0 \\end{matrix} \\right] =\n \\left[ \\begin{matrix} 3 \\\\ 2 \\\\ 0 \\end{matrix} \\right]\\end{align}\n\nFrom an information perspective, this projection has taken the point\n$x, y, z$ and removed the information about how far in the $z$\ndirection our point was located; all we know now is its position in the\n$x, y$ plane. Moreover, applying our projection matrix to *any point*\nin $x, y, z$ space will reduce it to a corresponding point on the\n$x, y$ plane. The term for this is a *subspace*: the projection matrix\nprojects points in the original space into a *subspace* of lower dimension\nthan the original. The reason our subspace is the $x,y$ plane (instead\nof, say, the $y,z$ plane) is a direct result of the particular values\nin our projection matrix.\n\n\n### Example: projection as noise reduction\n\nAnother way to describe this \"loss of information\" or \"projection into a\nsubspace\" is to say that projection reduces the rank (or \"degrees of\nfreedom\") of the measurement \u2014 here, from 3 dimensions down to 2. On the\nother hand, if you know that measurement component in the $z$ direction\nis just noise due to your measurement method, and all you care about are the\n$x$ and $y$ components, then projecting your 3-dimensional\nmeasurement into the $x, y$ plane could be seen as a form of noise\nreduction.\n\nOf course, it would be very lucky indeed if all the measurement noise were\nconcentrated in the $z$ direction; you could just discard the $z$\ncomponent without bothering to construct a projection matrix or do the matrix\nmultiplication. Suppose instead that in order to take that measurement you\nhad to pull a trigger on a measurement device, and the act of pulling the\ntrigger causes the device to move a little. If you measure how\ntrigger-pulling affects measurement device position, you could then \"correct\"\nyour real measurements to \"project out\" the effect of the trigger pulling.\nHere we'll suppose that the average effect of the trigger is to move the\nmeasurement device by $(3, -1, 1)$:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "trigger_effect = np.array([[3, -1, 1]]).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Knowing that, we can compute a plane that is orthogonal to the effect of the\ntrigger (using the fact that a plane through the origin has equation\n$Ax + By + Cz = 0$ given a normal vector $(A, B, C)$), and\nproject our real measurements onto that plane.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# compute the plane orthogonal to trigger_effect\nx, y = np.meshgrid(np.linspace(-1, 5, 61), np.linspace(-1, 5, 61))\nA, B, C = trigger_effect\nz = (-A * x - B * y) / C\n# cut off the plane below z=0 (just to make the plot nicer)\nmask = np.where(z >= 0)\nx = x[mask]\ny = y[mask]\nz = z[mask]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Computing the projection matrix from the trigger_effect vector is done\nusing [singular value decomposition](svd_) (SVD); interested readers may\nconsult the internet or a linear algebra textbook for details on this method.\nWith the projection matrix in place, we can project our original vector\n$(3, 2, 5)$ to remove the effect of the trigger, and then plot it:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# compute the projection matrix\nU, S, V = svd(trigger_effect, full_matrices=False)\ntrigger_projection_matrix = np.eye(3) - U @ U.T\n\n# project the vector onto the orthogonal plane\nprojected_point = trigger_projection_matrix @ point\nprojected_vector = trigger_projection_matrix @ vector\n\n# plot the trigger effect and its orthogonal plane\nax = setup_3d_axes()\nax.plot_trisurf(x, y, z, color='C2', shade=False, alpha=0.25)\nax.quiver3D(*np.concatenate([origin, trigger_effect]).flatten(),\n arrow_length_ratio=0.1, color='C2', alpha=0.5)\n\n# plot the original vector\nax.plot(*vector, color='k')\nax.plot(*point, color='k', marker='o')\noffset = np.full((3, 1), 0.1)\nax.text(*(point + offset).flat, '({}, {}, {})'.format(*point.flat), color='k')\n\n# plot the projected vector\nax.plot(*projected_vector, color='C0')\nax.plot(*projected_point, color='C0', marker='o')\noffset = np.full((3, 1), -0.2)\nax.text(*(projected_point + offset).flat,\n '({}, {}, {})'.format(*np.round(projected_point.flat, 2)),\n color='C0', horizontalalignment='right')\n\n# add dashed arrow showing projection\narrow_coords = np.concatenate([point, projected_point - point]).flatten()\nax.quiver3D(*arrow_coords, length=0.96, arrow_length_ratio=0.1,\n color='C1', linewidth=1, linestyle='dashed')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Just as before, the projection matrix will map *any point* in $x, y, z$\nspace onto that plane, and once a point has been projected onto that plane,\napplying the projection again will have no effect. For that reason, it should\nbe clear that although the projected points vary in all three $x$,\n$y$, and $z$ directions, the set of projected points have only\ntwo *effective* dimensions (i.e., they are constrained to a plane).\n\n.. sidebar:: Terminology\n\n In MNE-Python, the matrix used to project a raw signal into a subspace is\n usually called a :term:projector or a *projection\n operator* \u2014 these terms are interchangeable with the term *projection\n matrix* used above.\n\nProjections of EEG or MEG signals work in very much the same way: the point\n$x, y, z$ corresponds to the value of each sensor at a single time\npoint, and the projection matrix varies depending on what aspects of the\nsignal (i.e., what kind of noise) you are trying to project out. The only\nreal difference is that instead of a single 3-dimensional point $(x, y,\nz)$ you're dealing with a time series of $N$-dimensional \"points\" (one\nat each sampling time), where $N$ is usually in the tens or hundreds\n(depending on how many sensors your EEG/MEG system has). Fortunately, because\nprojection is a matrix operation, it can be done very quickly even on signals\nwith hundreds of dimensions and tens of thousands of time points.\n\n\n\n## Signal-space projection (SSP)\n\nWe mentioned above that the projection matrix will vary depending on what\nkind of noise you are trying to project away. Signal-space projection (SSP)\n:footcite:UusitaloIlmoniemi1997 is a way of estimating what that projection\nmatrix should be, by\ncomparing measurements with and without the signal of interest. For example,\nyou can take additional \"empty room\" measurements that record activity at the\nsensors when no subject is present. By looking at the spatial pattern of\nactivity across MEG sensors in an empty room measurement, you can create one\nor more $N$-dimensional vector(s) giving the \"direction(s)\" of\nenvironmental noise in sensor space (analogous to the vector for \"effect of\nthe trigger\" in our example above). SSP is also often used for removing\nheartbeat and eye movement artifacts \u2014 in those cases, instead of empty room\nrecordings the direction of the noise is estimated by detecting the\nartifacts, extracting epochs around them, and averaging. See\ntut-artifact-ssp for examples.\n\nOnce you know the noise vectors, you can create a hyperplane that is\northogonal\nto them, and construct a projection matrix to project your experimental\nrecordings onto that hyperplane. In that way, the component of your\nmeasurements associated with environmental noise can be removed. Again, it\nshould be clear that the projection reduces the dimensionality of your data \u2014\nyou'll still have the same number of sensor signals, but they won't all be\n*linearly independent* \u2014 but typically there are tens or hundreds of sensors\nand the noise subspace that you are eliminating has only 3-5 dimensions, so\nthe loss of degrees of freedom is usually not problematic.\n\n\n## Projectors in MNE-Python\n\nIn our example data, SSP  has already been performed\nusing empty room recordings, but the :term:projectors  are\nstored alongside the raw data and have not been *applied* yet (or,\nsynonymously, the projectors are not *active* yet). Here we'll load\nthe sample data  and crop it to 60 seconds; you can\nsee the projectors in the output of :func:~mne.io.read_raw_fif below:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sample_data_folder = mne.datasets.sample.data_path()\nsample_data_raw_file = os.path.join(sample_data_folder, 'MEG', 'sample',\n 'sample_audvis_raw.fif')\nraw = mne.io.read_raw_fif(sample_data_raw_file)\nraw.crop(tmax=60).load_data()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In MNE-Python, the environmental noise vectors are computed using [principal\ncomponent analysis](pca_), usually abbreviated \"PCA\", which is why the SSP\nprojectors usually have names like \"PCA-v1\". (Incidentally, since the process\nof performing PCA uses [singular value decomposition](svd_) under the hood,\nit is also common to see phrases like \"projectors were computed using SVD\" in\npublished papers.) The projectors are stored in the projs field of\nraw.info:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(raw.info['projs'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "raw.info['projs'] is an ordinary Python :class:list of\n:class:~mne.Projection objects, so you can access individual projectors by\nindexing into it. The :class:~mne.Projection object itself is similar to a\nPython :class:dict, so you can use its .keys() method to see what\nfields it contains (normally you don't need to access its properties\ndirectly, but you can if necessary):\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "first_projector = raw.info['projs']\nprint(first_projector)\nprint(first_projector.keys())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The :class:~mne.io.Raw, :class:~mne.Epochs, and :class:~mne.Evoked\nobjects all have a boolean :attr:~mne.io.Raw.proj attribute that indicates\nwhether there are any unapplied / inactive projectors stored in the object.\nIn other words, the :attr:~mne.io.Raw.proj attribute is True if at\nleast one :term:projector is present and all of them are active. In\naddition, each individual projector also has a boolean active field:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(raw.proj)\nprint(first_projector['active'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Computing projectors\n\nIn MNE-Python, SSP vectors can be computed using general purpose functions\n:func:mne.compute_proj_raw, :func:mne.compute_proj_epochs, and\n:func:mne.compute_proj_evoked. The general assumption these functions make\nis that the data passed contains raw data, epochs or averages of the artifact\nyou want to repair via projection. In practice this typically involves\ncontinuous raw data of empty room recordings or averaged ECG or EOG\nartifacts. A second set of high-level convenience functions is provided to\ncompute projection vectors for typical use cases. This includes\n:func:mne.preprocessing.compute_proj_ecg and\n:func:mne.preprocessing.compute_proj_eog for computing the ECG and EOG\nrelated artifact components, respectively; see tut-artifact-ssp for\nexamples of these uses. For computing the EEG reference signal as a\nprojector, the function :func:mne.set_eeg_reference can be used; see\ntut-set-eeg-ref for more information.\n\n

#### Warning

It is best to compute projectors only on channels that will be\n used (e.g., excluding bad channels). This ensures that\n projection vectors will remain ortho-normalized and that they\n properly capture the activity of interest.

\n\n\n### Visualizing the effect of projectors\n\nYou can see the effect the projectors are having on the measured signal by\ncomparing plots with and without the projectors applied. By default,\nraw.plot() will apply the projectors in the background before plotting\n(without modifying the :class:~mne.io.Raw object); you can control this\nwith the boolean proj parameter as shown below, or you can turn them on\nand off interactively with the projectors interface, accessed via the\n:kbd:Proj button in the lower right corner of the plot window. Here we'll\nlook at just the magnetometers, and a 2-second sample from the beginning of\nthe file.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mags = raw.copy().crop(tmax=2).pick_types(meg='mag')\nfor proj in (False, True):\n with mne.viz.use_browser_backend('matplotlib'):\n fig = mags.plot(butterfly=True, proj=proj)\n fig.subplots_adjust(top=0.9)\n fig.suptitle('proj={}'.format(proj), size='xx-large', weight='bold')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Additional ways of visualizing projectors are covered in the tutorial\ntut-artifact-ssp.\n\n\n### Loading and saving projectors\n\nSSP can be used for other types of signal cleaning besides just reduction of\nenvironmental noise. You probably noticed two large deflections in the\nmagnetometer signals in the previous plot that were not removed by the\nempty-room projectors \u2014 those are artifacts of the subject's heartbeat. SSP\ncan be used to remove those artifacts as well. The sample data includes\nprojectors for heartbeat noise reduction that were saved in a separate file\nfrom the raw data, which can be loaded with the :func:mne.read_proj\nfunction:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ecg_proj_file = os.path.join(sample_data_folder, 'MEG', 'sample',\n 'sample_audvis_ecg-proj.fif')\necg_projs = mne.read_proj(ecg_proj_file)\nprint(ecg_projs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is a corresponding :func:mne.write_proj function that can be used to\nsave projectors to disk in .fif format:\n\npython3\nmne.write_proj('heartbeat-proj.fif', ecg_projs)\n\n

#### Note

By convention, MNE-Python expects projectors to be saved with a filename\n ending in -proj.fif (or -proj.fif.gz), and will issue a warning\n if you forgo this recommendation.

\n\n\n### Adding and removing projectors\n\nAbove, when we printed the ecg_projs list that we loaded from a file, it\nshowed two projectors for gradiometers (the first two, marked \"planar\"), two\nfor magnetometers (the middle two, marked \"axial\"), and two for EEG sensors\n(the last two, marked \"eeg\"). We can add them to the :class:~mne.io.Raw\nobject using the :meth:~mne.io.Raw.add_proj method:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "raw.add_proj(ecg_projs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To remove projectors, there is a corresponding method\n:meth:~mne.io.Raw.del_proj that will remove projectors based on their index\nwithin the raw.info['projs'] list. For the special case of replacing the\nexisting projectors with new ones, use\nraw.add_proj(ecg_projs, remove_existing=True).\n\nTo see how the ECG projectors affect the measured signal, we can once again\nplot the data with and without the projectors applied (though remember that\nthe :meth:~mne.io.Raw.plot method only *temporarily* applies the projectors\nfor visualization, and does not permanently change the underlying data).\nWe'll compare the mags variable we created above, which had only the\nempty room SSP projectors, to the data with both empty room and ECG\nprojectors:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mags_ecg = raw.copy().crop(tmax=2).pick_types(meg='mag')\nfor data, title in zip([mags, mags_ecg], ['Without', 'With']):\n with mne.viz.use_browser_backend('matplotlib'):\n fig = data.plot(butterfly=True, proj=True)\n fig.subplots_adjust(top=0.9)\n fig.suptitle('{} ECG projector'.format(title), size='xx-large',\n weight='bold')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### When are projectors \"applied\"?\n\nBy default, projectors are applied when creating :class:epoched\n data from :class:~mne.io.Raw data, though application of the\nprojectors can be *delayed* by passing proj=False to the\n:class:~mne.Epochs constructor. However, even when projectors have not been\napplied, the :meth:mne.Epochs.get_data method will return data *as if the\nprojectors had been applied* (though the :class:~mne.Epochs object will be\nunchanged). Additionally, projectors cannot be applied if the data are not\npreloaded . If the data are memory-mapped_ (i.e., not\npreloaded), you can check the _projector attribute to see whether any\nprojectors will be applied once the data is loaded in memory.\n\nFinally, when performing inverse imaging (i.e., with\n:func:mne.minimum_norm.apply_inverse), the projectors will be\nautomatically applied. It is also possible to apply projectors manually when\nworking with :class:~mne.io.Raw, :class:~mne.Epochs or\n:class:~mne.Evoked objects via the object's :meth:~mne.io.Raw.apply_proj\nmethod. For all instance types, you can always copy the contents of\n:samp:{}.info['projs'] into a separate :class:list variable,\nuse :samp:{}.del_proj({}) to remove\none or more projectors, and then add them back later with\n:samp:{}.add_proj({}) if desired.\n\n

#### Warning

Remember that once a projector is applied, it can't be un-applied, so\n during interactive / exploratory analysis it's a good idea to use the\n object's :meth:~mne.io.Raw.copy method before applying projectors.

\n\n\n### Best practices\n\nIn general, it is recommended to apply projectors when creating\n:class:~mne.Epochs from :class:~mne.io.Raw data. There are two reasons\nfor this recommendation:\n\n1. It is computationally cheaper to apply projectors to data *after* the\n data have been reducted to just the segments of interest (the epochs)\n\n2. If you are applying amplitude-based rejection criteria to epochs, it is\n preferable to reject based on the signal *after* projectors have been\n applied, because the projectors may reduce noise in some epochs to\n tolerable levels (thereby increasing the number of acceptable epochs and\n consequenty increasing statistical power in any later analyses).\n\n\n## References\n\n.. footbibliography::\n\n\n.. LINKS\n\n https://docs.python.org/3/tutorial/controlflow.html#tut-unpacking-arguments\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.8.0" } }, "nbformat": 4, "nbformat_minor": 0 }