from os import path as op import numpy as np import matplotlib.pyplot as plt import mne from mne.forward import make_forward_solution from mne.io.constants import FIFF from mne.transforms import invert_transform, apply_trans from mne.source_space import setup_volume_source_space, read_source_spaces from mne.bem import read_bem_solution from mne.minimum_norm import make_inverse_operator, estimate_snr # imports needed for modified apply_inverse from mne.minimum_norm.inverse import (_check_reference, _check_method, _check_ori, _check_ch_names, prepare_inverse_operator, _pick_channels_inverse_operator, _assemble_kernel, combine_xyz, _subject_from_inverse, _make_stc, logger) def apply_inverse_mod(evoked, inverse_operator, lambda2=1. / 9., method="dSPM", pick_ori=None, prepared=False, label=None, verbose=None): """Apply inverse operator to evoked data. Modified to return the uncombined inverse solution. Parameters ---------- evoked : Evoked object Evoked data. inverse_operator: instance of InverseOperator Inverse operator returned from `mne.read_inverse_operator`, `prepare_inverse_operator` or `make_inverse_operator`. lambda2 : float The regularization parameter. method : "MNE" | "dSPM" | "sLORETA" Use mininum norm, dSPM or sLORETA. pick_ori : None | "normal" | "vector" If "normal", rather than pooling the orientations by taking the norm, only the radial component is kept. This is only implemented when working with loose orientations. If "vector", no pooling of the orientations is done and the vector result will be returned in the form of a :class:`mne.VectorSourceEstimate` object. This is only implemented when working with loose orientations. prepared : bool If True, do not call `prepare_inverse_operator`. label : Label | None Restricts the source estimates to a given label. If None, source estimates will be computed for the entire source space. verbose : bool, str, int, or None If not None, override default verbose level (see :func:`mne.verbose` and :ref:`Logging documentation ` for more). Returns ------- stc : SourceEstimate | VectorSourceEstimate | VolSourceEstimate The source estimates sol_orig: Array Results of applying the inverse kernel to evoked data See Also -------- apply_inverse_raw : Apply inverse operator to raw object apply_inverse_epochs : Apply inverse operator to epochs object """ _check_reference(evoked) _check_method(method) _check_ori(pick_ori, inverse_operator['source_ori']) # # Set up the inverse according to the parameters # nave = evoked.nave _check_ch_names(inverse_operator, evoked.info) if not prepared: inv = prepare_inverse_operator(inverse_operator, nave, lambda2, method) else: inv = inverse_operator # # Pick the correct channels from the data # sel = _pick_channels_inverse_operator(evoked.ch_names, inv) logger.info('Picked %d channels from the data' % len(sel)) logger.info('Computing inverse...') K, noise_norm, vertno, source_nn = _assemble_kernel(inv, label, method, pick_ori) sol = np.dot(K, evoked.data[sel]) # apply imaging kernel sol_orig = sol is_free_ori = (inverse_operator['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI and pick_ori != 'normal') if is_free_ori and pick_ori != 'vector': logger.info('combining the current components...') sol = combine_xyz(sol) if noise_norm is not None: logger.info('(dSPM)...') if is_free_ori and pick_ori == 'vector': noise_norm = noise_norm.repeat(3, axis=0) sol *= noise_norm tstep = 1.0 / evoked.info['sfreq'] tmin = float(evoked.times[0]) subject = _subject_from_inverse(inverse_operator) stc = _make_stc(sol, vertno, tmin=tmin, tstep=tstep, subject=subject, vector=(pick_ori == 'vector'), source_nn=source_nn) logger.info('[done]') return stc, sol_orig data_path = mne.datasets.sample.data_path() subjects_dir = op.join(data_path, 'subjects') fname_ave = op.join(data_path, 'MEG', 'sample', 'sample_audvis-ave.fif') fname_cov = op.join(data_path, 'MEG', 'sample', 'sample_audvis-cov.fif') fname_bem = op.join(subjects_dir, 'sample', 'bem', 'sample-5120-bem-sol.fif') fname_trans = op.join(data_path, 'MEG', 'sample', 'sample_audvis_raw-trans.fif') fname_src = op.join(subjects_dir, 'sample', 'bem', 'sample-oct-6-src.fif') bem = read_bem_solution(fname_bem) cov = mne.read_cov(fname_cov) src_cs = read_source_spaces(fname_src) ############################################################################### # Let's localize the N100m (using MEG only) evoked = mne.read_evokeds(fname_ave, condition='Right Auditory', baseline=(None, 0)) evoked.pick_types(meg=True, eeg=False) # localizet the N100 to get the best position evoked_N100 = evoked.copy() evoked_N100.crop(0.07, 0.08) # Fit a dipole dip = mne.fit_dipole(evoked_N100, fname_cov, fname_bem, fname_trans)[0] # Plot the result in 3D brain with the MRI image. dip.plot_locations(fname_trans, 'sample', subjects_dir, mode='orthoview') # find time point with highest GOF to plot best_idx = np.argmax(dip.gof) best_time = dip.times[best_idx] print('Highest GOF %0.1f%% at t=%0.1f ms with confidence volume %0.1f cm^3' % (dip.gof[best_idx], best_time * 1000, dip.conf['vol'][best_idx] * 100 ** 3)) # do a dipole fit from this location allowing the orientations to change dip_region = mne.fit_dipole(evoked, fname_cov, fname_bem, fname_trans, pos=dip.pos[best_idx])[0] ############################################################################### # read in transfprosm head_mri_trans = mne.read_trans(fname_trans) mri_head_trans = invert_transform(head_mri_trans) # conform transforms assert(mri_head_trans['to'] == FIFF.FIFFV_COORD_HEAD) assert(mri_head_trans['from'] == FIFF.FIFFV_COORD_MRI) ############################################################################### # fit anatomical source locations pos = dict() pos['rr'] = np.array([apply_trans(head_mri_trans, dip.pos[best_idx])]) # pos['nn'] = np.array([apply_trans(head_mri_trans, dip.ori[best_idx], False)]) pos['nn'] = np.array([[0, 0, 1]]) anat_src = setup_volume_source_space(subject=None, pos=pos, bem=fname_bem, mindist=5) # make forward operator fwd_anat = make_forward_solution(evoked.info, trans=fname_trans, src=anat_src, bem=bem, meg=True, eeg=False, n_jobs=1) if fwd_anat['nsource'] != len(anat_src[0]['rr']): raise RuntimeError('Some Sources were excluded!!!!') fwd_anat = mne.pick_types_forward(fwd_anat, meg=True, eeg=False) # make inverse operator inv_op_anat = make_inverse_operator(evoked.info, fwd_anat, cov, loose=1) lambda2s = [0.111, 0.5, 1., 3.16, 11.1] sol_orig_mag = dict() stc_anat_mag = dict() for l2 in lambda2s: stc_anat_mag[l2], sol_orig_mag[l2] = \ apply_inverse_mod(evoked, inv_op_anat, l2, method='dSPM', pick_ori=None) # illustrate the overall problem. plt.figure() for l2 in lambda2s: plt.plot(stc_anat_mag[l2].times, stc_anat_mag[l2].data[0], label='Amplitude (snr=%0.2f - l^2=%0.2f)' % (1./np.sqrt(l2), l2)) plt.plot(dip_region.times, dip_region.amplitude/dip_region.amplitude.max() * 35, lw=2, label='Regional dipole Amp (Scaled)', color='k') plt.legend() # SNR = 3 is just Bad how come? print ('Inverse operator singular Values (%0.2f, %0.2f, %0.2f)' % tuple(inv_op_anat['sing'])) # plot each orientation fig, ax = plt.subplots(3, 1) for l2 in lambda2s: for i in [0, 1, 2]: ax[i].plot(stc_anat_mag[l2].times, sol_orig_mag[l2][i], label='snr=%0.2f (l^2=%0.2f)' % (np.sqrt(1/l2), l2)) ax[i].set_title('Orientation %d' % i) for i in [0, 1, 2]: ax[i].legend() # orientation 1 is clearly very dependent on regularization # from the singular vaules I would have expected it to be orientation 2? # Can we use the data to predict the SNR to use? # compute snrs snr_whitened, snr_est = estimate_snr(evoked, inv_op_anat) fig, ax = plt.subplots(2, 1) ax[0].plot(evoked.times, snr_whitened, lw=2, label='SNR (Whitened)', color='r') ax[0].set_title('SNR estimated from Whitened Data') ax[1].plot(evoked.times, snr_est, lw=2, label='SNR (inv sol)', color='r') ax[1].set_title('SNR estimated inverse solutions') colors = plt.cm.get_cmap('jet', len(lambda2s)) for a in ax: for i, l2 in enumerate(lambda2s): a.axhline(y=np.sqrt(1/l2), label='snr=%0.2f (l^2=%0.2f)' % (np.sqrt(1/l2), l2), color=colors(i)) a.legend() ax[0].set_ylim(-0.1, 2 * snr_whitened.max()) ax[1].set_ylim(-0.1, 2 * snr_est.max()) plt.show() # yes