{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# KIT phantom dataset tutorial\n\nHere we read KIT data obtained from a phantom with 49 dipoles sequentially activated\nwith 2-cycle 11 Hz sinusoidal bursts to verify source localization accuracy.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Authors: Eric Larson \n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import mne_bids\nimport numpy as np\n\nimport mne\n\ndata_path = mne.datasets.phantom_kit.data_path()\nactual_pos, actual_ori = mne.dipole.get_phantom_dipoles(\"oyama\")\nactual_pos, actual_ori = actual_pos[:49], actual_ori[:49] # only 49 of 50 dipoles\n\nbids_path = mne_bids.BIDSPath(\n root=data_path,\n subject=\"01\",\n task=\"phantom\",\n run=\"01\",\n datatype=\"meg\",\n)\n# ignore warning about misc units\nraw = mne_bids.read_raw_bids(bids_path).load_data()\n\n# Let's apply a little bit of preprocessing (temporal filtering and reference\n# regression)\npicks_artifact = [\"MISC 001\", \"MISC 002\", \"MISC 003\"]\npicks = np.r_[\n mne.pick_types(raw.info, meg=True),\n mne.pick_channels(raw.info[\"ch_names\"], picks_artifact),\n]\nraw.filter(None, 40, picks=picks)\nmne.preprocessing.regress_artifact(\n raw, picks=\"meg\", picks_artifact=picks_artifact, copy=False, proj=False\n)\nplot_scalings = dict(mag=5e-12) # large-amplitude sinusoids\nraw_plot_kwargs = dict(duration=15, n_channels=50, scalings=plot_scalings)\nevents, event_id = mne.events_from_annotations(raw)\nraw.plot(events=events, **raw_plot_kwargs)\nn_dip = len(event_id)\nassert n_dip == 49 # sanity check" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also look at the power spectral density to see the phantom oscillations at\n11 Hz plus the expected frequency-domain sinc-like oscillations due to the time-domain\nboxcar windowing of the 11 Hz sinusoid.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "spectrum = raw.copy().crop(0, 60).compute_psd(n_fft=10000)\nfig = spectrum.plot(amplitude=False)\nfig.axes[0].set_xlim(0, 50)\ndip_freq = 11.0\nfig.axes[0].axvline(dip_freq, color=\"r\", ls=\"--\", lw=2, zorder=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can figure out our epoching parameters and epoch the data and plot it.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "tmin, tmax = -0.08, 0.18\nepochs = mne.Epochs(raw, tmin=tmin, tmax=tmax, decim=10, picks=\"data\", preload=True)\ndel raw\nepochs.plot(scalings=plot_scalings)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can average the epochs for each dipole, get the activation at the peak time,\nand create an :class:`mne.EvokedArray` from the result.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "t_peak = 1.0 / dip_freq / 4.0\ndata = np.zeros((len(epochs.ch_names), n_dip))\nfor di in range(n_dip):\n data[:, [di]] = epochs[f\"dip{di + 1:02d}\"].average().crop(t_peak, t_peak).data\nevoked = mne.EvokedArray(data, epochs.info, tmin=0, comment=\"KIT phantom activations\")\nevoked.plot_joint()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's fit dipoles at each dipole's peak activation time.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "trans = mne.transforms.Transform(\"head\", \"mri\", np.eye(4))\nsphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.08)\ncov = mne.compute_covariance(epochs, tmax=0, method=\"empirical\")\ndip, residual = mne.fit_dipole(evoked, cov, sphere, n_jobs=None)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally let's look at the results.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(f\"Average amplitude: {np.mean(dip.amplitude) * 1e9:0.1f} nAm\")\nprint(f\"Average GOF: {np.mean(dip.gof):0.1f}%\")\ndiffs = 1000 * np.sqrt(np.sum((dip.pos - actual_pos) ** 2, axis=-1))\nprint(f\"Average loc error: {np.mean(diffs):0.1f} mm\")\nangles = np.rad2deg(np.arccos(np.abs(np.sum(dip.ori * actual_ori, axis=1))))\nprint(f\"Average ori error: {np.mean(angles):0.1f}\u00b0\")\n\nfig = mne.viz.plot_alignment(\n evoked.info,\n trans,\n bem=sphere,\n coord_frame=\"head\",\n meg=\"helmet\",\n show_axes=True,\n)\nfig = mne.viz.plot_dipole_locations(\n dipoles=dip, mode=\"arrow\", color=(0.2, 1.0, 0.5), fig=fig\n)\n\nactual_amp = np.ones(len(dip)) # misc amp to create Dipole instance\nactual_gof = np.ones(len(dip)) # misc GOF to create Dipole instance\ndip_true = mne.Dipole(dip.times, actual_pos, actual_amp, actual_ori, actual_gof)\nfig = mne.viz.plot_dipole_locations(\n dipoles=dip_true, mode=\"arrow\", color=(0.0, 0.0, 0.0), fig=fig\n)\n\nmne.viz.set_3d_view(figure=fig, azimuth=90, elevation=90, distance=0.5)" ] } ], "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 }