Python mne.make_forward_solution函数代码示例

def compute_fwd_sol(raw_info, trans_fname, src, bem, fwd_filename):
    import mne

    mne.make_forward_solution(raw_info, trans_fname, src, bem,
                              fwd_filename,
                              mindist=5.0, # ignore sources <= 0mm from inner skull
                              meg=True, eeg=False,
                              n_jobs=2,
                              overwrite=True)

    print '\n*** FWD file %s written!!!\n' % fwd_filename

def test_inverse_operator_channel_ordering():
    """Test MNE inverse computation is immune to channel reorderings
    """
    # These are with original ordering
    evoked = _get_evoked()
    noise_cov = read_cov(fname_cov)

    fwd_orig = make_forward_solution(evoked.info, fname_trans, src_fname,
                                     fname_bem, eeg=True, mindist=5.0)
    fwd_orig = convert_forward_solution(fwd_orig, surf_ori=True)
    inv_orig = make_inverse_operator(evoked.info, fwd_orig, noise_cov,
                                     loose=0.2, depth=0.8,
                                     limit_depth_chs=False)
    stc_1 = apply_inverse(evoked, inv_orig, lambda2, "dSPM")

    # Assume that a raw reordering applies to both evoked and noise_cov,
    # so we don't need to create those from scratch. Just reorder them,
    # then try to apply the original inverse operator
    new_order = np.arange(len(evoked.info['ch_names']))
    randomiser = np.random.RandomState(42)
    randomiser.shuffle(new_order)
    evoked.data = evoked.data[new_order]
    evoked.info['chs'] = [evoked.info['chs'][n] for n in new_order]
    evoked.info._update_redundant()
    evoked.info._check_consistency()

    cov_ch_reorder = [c for c in evoked.info['ch_names']
                      if (c in noise_cov.ch_names)]

    new_order_cov = [noise_cov.ch_names.index(name) for name in cov_ch_reorder]
    noise_cov['data'] = noise_cov.data[np.ix_(new_order_cov, new_order_cov)]
    noise_cov['names'] = [noise_cov['names'][idx] for idx in new_order_cov]

    fwd_reorder = make_forward_solution(evoked.info, fname_trans, src_fname,
                                        fname_bem, eeg=True, mindist=5.0)
    fwd_reorder = convert_forward_solution(fwd_reorder, surf_ori=True)
    inv_reorder = make_inverse_operator(evoked.info, fwd_reorder, noise_cov,
                                        loose=0.2, depth=0.8,
                                        limit_depth_chs=False)

    stc_2 = apply_inverse(evoked, inv_reorder, lambda2, "dSPM")

    assert_equal(stc_1.subject, stc_2.subject)
    assert_array_equal(stc_1.times, stc_2.times)
    assert_allclose(stc_1.data, stc_2.data, rtol=1e-5, atol=1e-5)
    assert_true(inv_orig['units'] == inv_reorder['units'])

    # Reload with original ordering & apply reordered inverse
    evoked = _get_evoked()
    noise_cov = read_cov(fname_cov)

    stc_3 = apply_inverse(evoked, inv_reorder, lambda2, "dSPM")
    assert_allclose(stc_1.data, stc_3.data, rtol=1e-5, atol=1e-5)

def test_make_forward_solution():
    """Test making M-EEG forward solution from python."""
    fwd_py = make_forward_solution(fname_raw, fname_trans, fname_src,
                                   fname_bem, mindist=5.)
    assert (isinstance(fwd_py, Forward))
    fwd = read_forward_solution(fname_meeg)
    assert (isinstance(fwd, Forward))
    _compare_forwards(fwd, fwd_py, 366, 1494, meg_rtol=1e-3)
    # Homogeneous model
    with pytest.raises(RuntimeError, match='homogeneous.*1-layer.*EEG'):
        make_forward_solution(fname_raw, fname_trans, fname_src,
                              fname_bem_meg)

def test_simulate_round_trip():
    """Test simulate_raw round trip calculations."""
    # Check a diagonal round-trip
    raw, src, stc, trans, sphere = _get_data()
    raw.pick_types(meg=True, stim=True)
    bem = read_bem_solution(bem_1_fname)
    old_bem = bem.copy()
    old_src = src.copy()
    old_trans = trans.copy()
    fwd = make_forward_solution(raw.info, trans, src, bem)
    # no omissions
    assert (sum(len(s['vertno']) for s in src) ==
            sum(len(s['vertno']) for s in fwd['src']) ==
            36)
    # make sure things were not modified
    assert (old_bem['surfs'][0]['coord_frame'] ==
            bem['surfs'][0]['coord_frame'])
    assert trans == old_trans
    _compare_source_spaces(src, old_src)
    data = np.eye(fwd['nsource'])
    raw.crop(0, (len(data) - 1) / raw.info['sfreq'])
    stc = SourceEstimate(data, [s['vertno'] for s in fwd['src']],
                         0, 1. / raw.info['sfreq'])
    for use_cps in (False, True):
        this_raw = simulate_raw(raw, stc, trans, src, bem, cov=None,
                                use_cps=use_cps)
        this_raw.pick_types(meg=True, eeg=True)
        assert (old_bem['surfs'][0]['coord_frame'] ==
                bem['surfs'][0]['coord_frame'])
        assert trans == old_trans
        _compare_source_spaces(src, old_src)
        this_fwd = convert_forward_solution(fwd, force_fixed=True,
                                            use_cps=use_cps)
        assert_allclose(this_raw[:][0], this_fwd['sol']['data'],
                        atol=1e-12, rtol=1e-6)

def test_gamma_map_vol_sphere():
    """Gamma MAP with a sphere forward and volumic source space"""
    evoked = read_evokeds(fname_evoked, condition=0, baseline=(None, 0),
                          proj=False)
    evoked.resample(50, npad=100)
    evoked.crop(tmin=0.1, tmax=0.16)  # crop to window around peak

    cov = read_cov(fname_cov)
    cov = regularize(cov, evoked.info)

    info = evoked.info
    sphere = mne.make_sphere_model(r0=(0., 0., 0.), head_radius=0.080)
    src = mne.setup_volume_source_space(subject=None, pos=15., mri=None,
                                        sphere=(0.0, 0.0, 0.0, 80.0),
                                        bem=None, mindist=5.0,
                                        exclude=2.0)
    fwd = mne.make_forward_solution(info, trans=None, src=src, bem=sphere,
                                    eeg=False, meg=True)

    alpha = 0.5
    assert_raises(ValueError, gamma_map, evoked, fwd, cov, alpha,
                  loose=0, return_residual=False)

    assert_raises(ValueError, gamma_map, evoked, fwd, cov, alpha,
                  loose=0.2, return_residual=False)

    stc = gamma_map(evoked, fwd, cov, alpha, tol=1e-4,
                    xyz_same_gamma=False, update_mode=2,
                    return_residual=False)

    assert_array_almost_equal(stc.times, evoked.times, 5)

def run_forward(subject_id):
    subject = "sub%03d" % subject_id
    print("processing subject: %s" % subject)
    data_path = op.join(meg_dir, subject)

    fname_ave = op.join(data_path, '%s-ave.fif' % subject)
    fname_fwd = op.join(data_path, '%s-meg-%s-fwd.fif' % (subject, spacing))
    fname_trans = op.join(study_path, 'ds117', subject, 'MEG', '%s-trans.fif' % subject)

    src = mne.setup_source_space(subject, spacing=spacing,
                                 subjects_dir=subjects_dir, overwrite=True,
                                 n_jobs=1, add_dist=False)

    src_fname = op.join(subjects_dir, subject, '%s-src.fif' % spacing)
    mne.write_source_spaces(src_fname, src)

    bem_model = mne.make_bem_model(subject, ico=4, subjects_dir=subjects_dir,
                                   conductivity=(0.3,))
    bem = mne.make_bem_solution(bem_model)
    info = mne.read_evokeds(fname_ave, condition=0).info
    fwd = mne.make_forward_solution(info, trans=fname_trans, src=src, bem=bem,
                                    fname=None, meg=True, eeg=False,
                                    mindist=mindist, n_jobs=1, overwrite=True)
    fwd = mne.convert_forward_solution(fwd, surf_ori=True)
    mne.write_forward_solution(fname_fwd, fwd, overwrite=True)

def test_make_forward_solution_sphere():
    """Test making a forward solution with a sphere model"""
    temp_dir = _TempDir()
    fname_src_small = op.join(temp_dir, "sample-oct-2-src.fif")
    src = setup_source_space("sample", fname_src_small, "oct2", subjects_dir=subjects_dir, add_dist=False)
    out_name = op.join(temp_dir, "tmp-fwd.fif")
    run_subprocess(
        [
            "mne_forward_solution",
            "--meg",
            "--eeg",
            "--meas",
            fname_raw,
            "--src",
            fname_src_small,
            "--mri",
            fname_trans,
            "--fwd",
            out_name,
        ]
    )
    fwd = read_forward_solution(out_name)
    sphere = make_sphere_model(verbose=True)
    fwd_py = make_forward_solution(fname_raw, fname_trans, src, sphere, meg=True, eeg=True, verbose=True)
    _compare_forwards(fwd, fwd_py, 366, 108, meg_rtol=5e-1, meg_atol=1e-6, eeg_rtol=5e-1, eeg_atol=5e-1)
    # Since the above is pretty lax, let's check a different way
    for meg, eeg in zip([True, False], [False, True]):
        fwd_ = pick_types_forward(fwd, meg=meg, eeg=eeg)
        fwd_py_ = pick_types_forward(fwd, meg=meg, eeg=eeg)
        assert_allclose(np.corrcoef(fwd_["sol"]["data"].ravel(), fwd_py_["sol"]["data"].ravel())[0, 1], 1.0, rtol=1e-3)

def test_apply_inverse_sphere():
    """Test applying an inverse with a sphere model (rank-deficient)."""
    evoked = _get_evoked()
    evoked.pick_channels(evoked.ch_names[:306:8])
    evoked.info['projs'] = []
    cov = make_ad_hoc_cov(evoked.info)
    sphere = make_sphere_model('auto', 'auto', evoked.info)
    fwd = read_forward_solution(fname_fwd)
    vertices = [fwd['src'][0]['vertno'][::5],
                fwd['src'][1]['vertno'][::5]]
    stc = SourceEstimate(np.zeros((sum(len(v) for v in vertices), 1)),
                         vertices, 0., 1.)
    fwd = restrict_forward_to_stc(fwd, stc)
    fwd = make_forward_solution(evoked.info, fwd['mri_head_t'], fwd['src'],
                                sphere, mindist=5.)
    evoked = EvokedArray(fwd['sol']['data'].copy(), evoked.info)
    assert fwd['sol']['nrow'] == 39
    assert fwd['nsource'] == 101
    assert fwd['sol']['ncol'] == 303
    tempdir = _TempDir()
    temp_fname = op.join(tempdir, 'temp-inv.fif')
    inv = make_inverse_operator(evoked.info, fwd, cov, loose=1.)
    # This forces everything to be float32
    write_inverse_operator(temp_fname, inv)
    inv = read_inverse_operator(temp_fname)
    stc = apply_inverse(evoked, inv, method='eLORETA',
                        method_params=dict(eps=1e-2))
    # assert zero localization bias
    assert_array_equal(np.argmax(stc.data, axis=0),
                       np.repeat(np.arange(101), 3))

def test_make_forward_solution_sphere():
    """Test making a forward solution with a sphere model."""
    temp_dir = _TempDir()
    fname_src_small = op.join(temp_dir, 'sample-oct-2-src.fif')
    src = setup_source_space('sample', 'oct2', subjects_dir=subjects_dir,
                             add_dist=False)
    write_source_spaces(fname_src_small, src)  # to enable working with MNE-C
    out_name = op.join(temp_dir, 'tmp-fwd.fif')
    run_subprocess(['mne_forward_solution', '--meg', '--eeg',
                    '--meas', fname_raw, '--src', fname_src_small,
                    '--mri', fname_trans, '--fwd', out_name])
    fwd = read_forward_solution(out_name)
    sphere = make_sphere_model(verbose=True)
    fwd_py = make_forward_solution(fname_raw, fname_trans, src, sphere,
                                   meg=True, eeg=True, verbose=True)
    _compare_forwards(fwd, fwd_py, 366, 108,
                      meg_rtol=5e-1, meg_atol=1e-6,
                      eeg_rtol=5e-1, eeg_atol=5e-1)
    # Since the above is pretty lax, let's check a different way
    for meg, eeg in zip([True, False], [False, True]):
        fwd_ = pick_types_forward(fwd, meg=meg, eeg=eeg)
        fwd_py_ = pick_types_forward(fwd, meg=meg, eeg=eeg)
        assert_allclose(np.corrcoef(fwd_['sol']['data'].ravel(),
                                    fwd_py_['sol']['data'].ravel())[0, 1],
                        1.0, rtol=1e-3)

    def compute_forward_and_inverse_solutions(self, orientation = 'fixed'):
        """docstring for compute_forward_solution"""

        info = self.grand_average_evoked.info
        trans = mne.read_trans(op.join(self.processed_files, '%s-trans.fif' %self.subject))
        src = glob.glob(op.join(self.subjects_dir, self.subject, 'bem', '*-ico-4-src.fif'))[0]
        bem = glob.glob(op.join(self.subjects_dir, self.subject, 'bem', '*-bem-sol.fif'))[0]
        fname = op.join(self.processed_files, '%s_forward.fif' %self.subject)

        # check if fwd exists, if not, make it
        if not op.exists(fname):
            fwd = mne.make_forward_solution(info = info, trans = trans, src = src,
                                            bem = bem, fname = fname, meg = True, eeg = False,
                                            overwrite = True, ignore_ref = True)

            self.add_preprocessing_notes("Forward solution generated and saved to %s" %fname)

        if orientation == 'fixed':
            force_fixed = True
        else:
            force_fixed = False

        fwd = mne.read_forward_solution(fname,force_fixed=force_fixed)

        self.forward_solution = fwd

        inv = mne.minimum_norm.make_inverse_operator(info, self.forward_solution, self.cov_reg, loose = None, depth = None, fixed = force_fixed)
        self.inverse_solution = inv
        mne.minimum_norm.write_inverse_operator(op.join(self.processed_files, '%s_inv.fif' %self.subject), self.inverse_solution)

        self.add_preprocessing_notes("Inverse solution generated and saved to %s" %op.join(self.processed_files, '%s_inv.fif' %self.subject))
        return fwd, inv

def test_volume_labels_morph(tmpdir):
    """Test generating a source space from volume label."""
    # see gh-5224
    evoked = mne.read_evokeds(fname_evoked)[0].crop(0, 0)
    evoked.pick_channels(evoked.ch_names[:306:8])
    evoked.info.normalize_proj()
    n_ch = len(evoked.ch_names)
    aseg_fname = op.join(subjects_dir, 'sample', 'mri', 'aseg.mgz')
    label_names = get_volume_labels_from_aseg(aseg_fname)
    src = setup_volume_source_space(
        'sample', subjects_dir=subjects_dir, volume_label=label_names[:2],
        mri=aseg_fname)
    assert len(src) == 2
    assert src.kind == 'volume'
    n_src = sum(s['nuse'] for s in src)
    sphere = make_sphere_model('auto', 'auto', evoked.info)
    fwd = make_forward_solution(evoked.info, fname_trans, src, sphere)
    assert fwd['sol']['data'].shape == (n_ch, n_src * 3)
    inv = make_inverse_operator(evoked.info, fwd, make_ad_hoc_cov(evoked.info),
                                loose=1.)
    stc = apply_inverse(evoked, inv)
    assert stc.data.shape == (n_src, 1)
    img = stc.as_volume(src, mri_resolution=True)
    n_on = np.array(img.dataobj).astype(bool).sum()
    assert n_on == 291  # was 291 on `master` before gh-5590
    img = stc.as_volume(src, mri_resolution=False)
    n_on = np.array(img.dataobj).astype(bool).sum()
    assert n_on == 44  # was 20 on `master` before gh-5590

def apply_inverse_ave(fnevo, min_subject='fsaverage'):
    
    from mne import make_forward_solution
    from mne.minimum_norm import write_inverse_operator
    fnlist = get_files_from_list(fnevo)
    # loop across all filenames
    for fname in fnlist:
        fn_path = os.path.split(fname)[0]
        name = os.path.basename(fname)
        #fn_inv = fname[:fname.rfind('-ave.fif')] + ',ave-inv.fif' 
        subject = name.split('_')[0]
        fn_inv = fn_path + '/%s_fibp1-45,ave-inv.fif' %subject
        subject_path = subjects_dir + '/%s' %subject
        #min_dir = subjects_dir + '/%s' %min_subject
        fn_trans = fn_path + '/%s-trans.fif' % subject
        #fn_cov = fn_path + '/%s_empty,nr,fibp1-45-cov.fif' % subject
        fn_cov = fn_path + '/%s_empty,fibp1-45-cov.fif' %subject
        fn_src = subject_path + '/bem/%s-oct-6-src.fif' % subject
        fn_bem = subject_path + '/bem/%s-5120-5120-5120-bem-sol.fif' % subject
        [evoked] = mne.read_evokeds(fname)
        evoked.pick_types(meg=True, ref_meg=False)
        noise_cov = mne.read_cov(fn_cov)
        #noise_cov = mne.cov.regularize(noise_cov, evoked.info,
         #                               mag=0.05, grad=0.05, proj=True)
        fwd = make_forward_solution(evoked.info, fn_trans, fn_src, fn_bem)
        fwd['surf_ori'] = True
        inv = mne.minimum_norm.make_inverse_operator(evoked.info, fwd, noise_cov, loose=0.2,
                                                     depth=0.8, limit_depth_chs=False)
        write_inverse_operator(fn_inv, inv)

def apply_inverse(fnevo, method='dSPM', snr=3.0, event='LLst', 
                  baseline=False, btmin=-0.3, btmax=-0.1, min_subject='fsaverage'):
    '''  
        Parameter
        ---------
        fnevo: string or list
            The evoked file with ECG, EOG and environmental noise free.
        method: inverse method, 'MNE' or 'dSPM'
        event: string
            The event name related with epochs.
        min_subject: string
            The subject name as the common brain.
        snr: signal to noise ratio for inverse solution. 
    '''
    #Get the default subjects_dir
    from mne.minimum_norm import apply_inverse
    fnlist = get_files_from_list(fnevo)
    # loop across all filenames
    for fname in fnlist:
        fn_path = os.path.split(fname)[0]
        name = os.path.basename(fname)
        stc_name = name[:name.rfind('-ave.fif')] 
        subject = name.split('_')[0]
        subject_path = subjects_dir + '/%s' %subject
        min_dir = subjects_dir + '/%s' %min_subject
        fn_trans = fn_path + '/%s-trans.fif' % subject
        fn_cov = fn_path + '/%s_empty,nr-cov.fif' % subject
        fn_src = subject_path + '/bem/%s-ico-5-src.fif' % subject
        fn_bem = subject_path + '/bem/%s-5120-5120-5120-bem-sol.fif' % subject
        snr = snr
        lambda2 = 1.0 / snr ** 2 
        #noise_cov = mne.read_cov(fn_cov)
        [evoked] = mne.read_evokeds(fname)
        noise_cov = mne.read_cov(fn_cov)
        # this path used for ROI definition
        stc_path = min_dir + '/%s_ROIs/%s' %(method,subject)
        #fn_cov = meg_path + '/%s_empty,fibp1-45,nr-cov.fif' % subject
        set_directory(stc_path)
        noise_cov = mne.cov.regularize(noise_cov, evoked.info,
                                        mag=0.05, grad=0.05, proj=True)
        fwd_ev = mne.make_forward_solution(evoked.info, trans=fn_trans,
                                            src=fn_src, bem=fn_bem,
                                            fname=None, meg=True, eeg=False,
                                            mindist=5.0, n_jobs=2,
                                            overwrite=True)
        fwd_ev = mne.convert_forward_solution(fwd_ev, surf_ori=True)
        forward_meg_ev = mne.pick_types_forward(fwd_ev, meg=True, eeg=False)
        inverse_operator_ev = mne.minimum_norm.make_inverse_operator(
            evoked.info, forward_meg_ev, noise_cov,
            loose=0.2, depth=0.8)
        # Compute inverse solution
        stc = apply_inverse(evoked, inverse_operator_ev, lambda2, method,
                            pick_ori=None)
        # Morph STC
        stc_morph = mne.morph_data(subject, min_subject, stc, grade=5, smooth=5)
        stc_morph.save(stc_path + '/%s' % (stc_name), ftype='stc')
        if baseline == True:
            stc_base = stc_morph.crop(btmin, btmax)
            stc_base.save(stc_path + '/%s_%s_baseline' % (subject, event), ftype='stc')

def test_simulate_round_trip(raw_data):
    """Test simulate_raw round trip calculations."""
    # Check a diagonal round-trip
    raw, src, stc, trans, sphere = raw_data
    raw.pick_types(meg=True, stim=True)
    bem = read_bem_solution(bem_1_fname)
    old_bem = bem.copy()
    old_src = src.copy()
    old_trans = trans.copy()
    fwd = make_forward_solution(raw.info, trans, src, bem)
    # no omissions
    assert (sum(len(s['vertno']) for s in src) ==
            sum(len(s['vertno']) for s in fwd['src']) ==
            36)
    # make sure things were not modified
    assert (old_bem['surfs'][0]['coord_frame'] ==
            bem['surfs'][0]['coord_frame'])
    assert trans == old_trans
    _compare_source_spaces(src, old_src)
    data = np.eye(fwd['nsource'])
    raw.crop(0, (len(data) - 1) / raw.info['sfreq'])
    stc = SourceEstimate(data, [s['vertno'] for s in fwd['src']],
                         0, 1. / raw.info['sfreq'])
    for use_fwd in (None, fwd):
        if use_fwd is None:
            use_trans, use_src, use_bem = trans, src, bem
        else:
            use_trans = use_src = use_bem = None
        with pytest.deprecated_call():
            this_raw = simulate_raw(raw, stc, use_trans, use_src, use_bem,
                                    cov=None, forward=use_fwd)
        this_raw.pick_types(meg=True, eeg=True)
        assert (old_bem['surfs'][0]['coord_frame'] ==
                bem['surfs'][0]['coord_frame'])
        assert trans == old_trans
        _compare_source_spaces(src, old_src)
        this_fwd = convert_forward_solution(fwd, force_fixed=True)
        assert_allclose(this_raw[:][0], this_fwd['sol']['data'],
                        atol=1e-12, rtol=1e-6)
    with pytest.raises(ValueError, match='If forward is not None then'):
        simulate_raw(raw.info, stc, trans, src, bem, forward=fwd)
    # Not iterable
    with pytest.raises(TypeError, match='SourceEstimate, tuple, or iterable'):
        simulate_raw(raw.info, 0., trans, src, bem, None)
    # STC with a source that `src` does not have
    assert 0 not in src[0]['vertno']
    vertices = [[0, fwd['src'][0]['vertno'][0]], []]
    stc_bad = SourceEstimate(data[:2], vertices, 0, 1. / raw.info['sfreq'])
    with pytest.warns(RuntimeWarning,
                      match='1 of 2 SourceEstimate vertices'):
        simulate_raw(raw.info, stc_bad, trans, src, bem)
    assert 0 not in fwd['src'][0]['vertno']
    with pytest.warns(RuntimeWarning,
                      match='1 of 2 SourceEstimate vertices'):
        simulate_raw(raw.info, stc_bad, None, None, None, forward=fwd)
    # dev_head_t mismatch
    fwd['info']['dev_head_t']['trans'][0, 0] = 1.
    with pytest.raises(ValueError, match='dev_head_t.*does not match'):
        simulate_raw(raw.info, stc, None, None, None, forward=fwd)

def test_make_forward_solution():
    """Test making M-EEG forward solution from python."""
    fwd_py = make_forward_solution(fname_raw, fname_trans, fname_src,
                                   fname_bem, mindist=5.0, eeg=True, meg=True)
    assert (isinstance(fwd_py, Forward))
    fwd = read_forward_solution(fname_meeg)
    assert (isinstance(fwd, Forward))
    _compare_forwards(fwd, fwd_py, 366, 1494, meg_rtol=1e-3)