Python transforms.apply_trans函数代码示例

def test_set_dig_montage():
    """Test applying DigMontage to inst."""
    # Extensive testing of applying `dig` to info is done in test_meas_info
    # with `test_make_dig_points`.
    names = ["nasion", "lpa", "rpa", "1", "2", "3", "4", "5"]
    hsp_points = _read_dig_points(hsp)
    elp_points = _read_dig_points(elp)
    nasion, lpa, rpa = elp_points[:3]
    nm_trans = get_ras_to_neuromag_trans(nasion, lpa, rpa)
    elp_points = apply_trans(nm_trans, elp_points)
    nasion, lpa, rpa = elp_points[:3]
    hsp_points = apply_trans(nm_trans, hsp_points)

    montage = read_dig_montage(hsp, hpi, elp, names, transform=True, dev_head_t=True)
    temp_dir = _TempDir()
    fname_temp = op.join(temp_dir, "test.fif")
    montage.save(fname_temp)
    montage_read = read_dig_montage(fif=fname_temp)
    for use_mon in (montage, montage_read):
        info = create_info(["Test Ch"], 1e3, ["eeg"])
        with warnings.catch_warnings(record=True) as w:  # test ch pos not set
            _set_montage(info, use_mon)
        assert_true(all("not set" in str(ww.message) for ww in w))
        hs = np.array([p["r"] for i, p in enumerate(info["dig"]) if p["kind"] == FIFF.FIFFV_POINT_EXTRA])
        nasion_dig = np.array(
            [
                p["r"]
                for p in info["dig"]
                if all([p["ident"] == FIFF.FIFFV_POINT_NASION, p["kind"] == FIFF.FIFFV_POINT_CARDINAL])
            ]
        )
        lpa_dig = np.array(
            [
                p["r"]
                for p in info["dig"]
                if all([p["ident"] == FIFF.FIFFV_POINT_LPA, p["kind"] == FIFF.FIFFV_POINT_CARDINAL])
            ]
        )
        rpa_dig = np.array(
            [
                p["r"]
                for p in info["dig"]
                if all([p["ident"] == FIFF.FIFFV_POINT_RPA, p["kind"] == FIFF.FIFFV_POINT_CARDINAL])
            ]
        )
        hpi_dig = np.array([p["r"] for p in info["dig"] if p["kind"] == FIFF.FIFFV_POINT_HPI])
        assert_allclose(hs, hsp_points, atol=1e-7)
        assert_allclose(nasion_dig.ravel(), nasion, atol=1e-7)
        assert_allclose(lpa_dig.ravel(), lpa, atol=1e-7)
        assert_allclose(rpa_dig.ravel(), rpa, atol=1e-7)
        assert_allclose(hpi_dig, elp_points[3:], atol=1e-7)

def test_fit_matched_points():
    """Test fit_matched_points: fitting two matching sets of points"""
    tgt_pts = np.random.uniform(size=(6, 3))

    # rotation only
    trans = rotation(2, 6, 3)
    src_pts = apply_trans(trans, tgt_pts)
    trans_est = fit_matched_points(src_pts, tgt_pts, translate=False, out="trans")
    est_pts = apply_trans(trans_est, src_pts)
    assert_array_almost_equal(tgt_pts, est_pts, 2, "fit_matched_points with " "rotation")

    # rotation & scaling
    trans = np.dot(rotation(2, 6, 3), scaling(0.5, 0.5, 0.5))
    src_pts = apply_trans(trans, tgt_pts)
    trans_est = fit_matched_points(src_pts, tgt_pts, translate=False, scale=1, out="trans")
    est_pts = apply_trans(trans_est, src_pts)
    assert_array_almost_equal(tgt_pts, est_pts, 2, "fit_matched_points with " "rotation and scaling.")

    # rotation & translation
    trans = np.dot(translation(2, -6, 3), rotation(2, 6, 3))
    src_pts = apply_trans(trans, tgt_pts)
    trans_est = fit_matched_points(src_pts, tgt_pts, out="trans")
    est_pts = apply_trans(trans_est, src_pts)
    assert_array_almost_equal(tgt_pts, est_pts, 2, "fit_matched_points with " "rotation and translation.")

    # rotation & translation & scaling
    trans = reduce(np.dot, (translation(2, -6, 3), rotation(1.5, 0.3, 1.4), scaling(0.5, 0.5, 0.5)))
    src_pts = apply_trans(trans, tgt_pts)
    trans_est = fit_matched_points(src_pts, tgt_pts, scale=1, out="trans")
    est_pts = apply_trans(trans_est, src_pts)
    assert_array_almost_equal(tgt_pts, est_pts, 2, "fit_matched_points with " "rotation, translation and scaling.")

    # test exceeding tolerance
    tgt_pts[0, :] += 20
    assert_raises(RuntimeError, fit_matched_points, tgt_pts, src_pts, tol=10)

def plot_coregistration(subject, subjects_dir, hcp_path, recordings_path,
                        info_from=(('data_type', 'rest'), ('run_index', 0)),
                        view_init=(('azim', 0), ('elev', 0))):
    """A diagnostic plot to show the HCP coregistration

    Parameters
    ----------
    subject : str
        The subject
    subjects_dir : str
        The path corresponding to MNE/freesurfer SUBJECTS_DIR (to be created)
    hcp_path : str
        The path where the HCP files can be found.
    recordings_path : str
        The path to converted data (including the head<->device transform).
    info_from : tuple of tuples | dict
        The reader info concerning the data from which sensor positions
        should be read.
        Must not be empty room as sensor positions are in head
        coordinates for 4D systems, hence not available in that case.
        Note that differences between the sensor positions across runs
        are smaller than 12 digits, hence negligible.
    view_init : tuple of tuples | dict
        The initival view, defaults to azimuth and elevation of 0,
        a simple lateral view

    Returns
    -------
    fig : matplotlib.figure.Figure
        The figure object.
    """
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D  #  noqa

    if isinstance(info_from, tuple):
        info_from = dict(info_from)
    if isinstance(view_init, tuple):
        view_init = dict(view_init)

    head_mri_t = read_trans(
        op.join(recordings_path, subject,
                '{}-head_mri-trans.fif'.format(subject)))

    info = read_info(subject=subject, hcp_path=hcp_path, **info_from)

    info = pick_info(info, _pick_data_channels(info, with_ref_meg=False))
    sens_pnts = np.array([c['loc'][:3] for c in info['chs']])
    sens_pnts = apply_trans(head_mri_t, sens_pnts)
    sens_pnts *= 1e3  # put in mm scale

    pnts, tris = read_surface(
        op.join(subjects_dir, subject, 'bem', 'inner_skull.surf'))

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(*sens_pnts.T, color='purple', marker='o')
    ax.scatter(*pnts.T, color='green', alpha=0.3)
    ax.view_init(**view_init)
    fig.tight_layout()
    return fig

#.........这里部分代码省略.........
    hcp_path : str
        The path where the HCP files can be found.
    outputs : {'label', 'mri', 'stats', 'surf', 'touch'}
        The outputs of the freesrufer pipeline shipped by HCP. Defaults to
        ('mri', 'surf'), the minimum needed to extract MNE-friendly anatomy
        files and data.
    """
    if hcp_path == op.curdir:
        hcp_path = op.realpath(hcp_path)
    if not op.isabs(subjects_dir):
        subjects_dir = op.realpath(subjects_dir)

    this_subjects_dir = op.join(subjects_dir, subject)
    if not op.isabs(recordings_path):
        recordings_path = op.realpath(recordings_path)

    this_recordings_path = op.join(recordings_path, subject)

    if not op.exists(this_recordings_path):
        os.makedirs(this_recordings_path)

    for output in outputs:
        if not op.exists(op.join(this_subjects_dir, output)):
            os.makedirs(op.join(this_subjects_dir, output))
        if output == 'mri':
            for suboutput in ['orig', 'transforms']:
                if not op.exists(
                        op.join(this_subjects_dir, output, suboutput)):
                    os.makedirs(op.join(this_subjects_dir, output, suboutput))

        files = get_file_paths(
            subject=subject, data_type='freesurfer', output=output,
            hcp_path=hcp_path)
        for source in files:
            match = [match for match in re.finditer(subject, source)][-1]
            split_path = source[:match.span()[1] + 1]
            target = op.join(this_subjects_dir, source.split(split_path)[-1])
            if (not op.isfile(target) and not op.islink(target) and
                    op.exists(source)):  # don't link if it's not there.
                if sys.platform != 'win32':
                    os.symlink(source, target)
                else:
                    shutil.copyfile(source, target)

    logger.info('reading extended structural processing ...')

    # Step 1 #################################################################
    # transform head models to expected coordinate system

    # make hcp trans
    transforms_fname = get_file_paths(
        subject=subject, data_type='meg_anatomy', output='transforms',
        hcp_path=hcp_path)
    transforms_fname = [k for k in transforms_fname if
                        k.endswith('transform.txt')][0]
    hcp_trans = _read_trans_hcp(fname=transforms_fname, convert_to_meter=False)

    # get RAS freesurfer trans
    c_ras_trans_fname = get_file_paths(
        subject=subject, data_type='freesurfer', output='mri',
        hcp_path=hcp_path)
    c_ras_trans_fname = [k for k in c_ras_trans_fname if
                         k.endswith('c_ras.mat')][0]
    logger.info('reading RAS freesurfer transform')
    # ceci n'est pas un .mat file ...

    with open(op.join(subjects_dir, c_ras_trans_fname)) as fid:
        ras_trans = np.array([
            r.split() for r in fid.read().split('\n') if r],
            dtype=np.float64)

    logger.info('Combining RAS transform and coregistration')
    ras_trans_m = linalg.inv(ras_trans)  # and the inversion

    logger.info('extracting head model')
    head_model_fname = get_file_paths(
        subject=subject, data_type='meg_anatomy', output='head_model',
        hcp_path=hcp_path)[0]
    pnts, faces = _get_head_model(head_model_fname=head_model_fname)

    logger.info('coregistring head model to MNE-HCP coordinates')
    pnts = apply_trans(ras_trans_m.dot(hcp_trans['bti2spm']), pnts)

    tri_fname = op.join(this_subjects_dir, 'bem', 'inner_skull.surf')
    if not op.exists(op.dirname(tri_fname)):
        os.makedirs(op.dirname(tri_fname))
    write_surface(tri_fname, pnts, faces)

    # Step 2 #################################################################
    # write corresponding device to MRI transform

    logger.info('extracting coregistration')
    # now convert to everything meter too here
    ras_trans_m[:3, 3] *= 1e-3
    bti2spm = hcp_trans['bti2spm']
    bti2spm[:3, 3] *= 1e-3
    head_mri_t = Transform(  # we're lying here for a good purpose
        'head', 'mri', np.dot(ras_trans_m, bti2spm))  # it should be 'ctf_head'
    write_trans(op.join(this_recordings_path,
                        '%s-head_mri-trans.fif' % subject), head_mri_t)

def _compute_depth(dip, fname_bem, fname_trans, subject, subjects_dir):
    """Compute dipole depth."""
    trans = _get_trans(fname_trans)[0]
    bem = read_bem_solution(fname_bem)
    surf = _bem_find_surface(bem, 'inner_skull')
    points = surf['rr']
    points = apply_trans(trans['trans'], points)
    depth = _compute_nearest(points, dip.pos, return_dists=True)[1][0]
    return np.ravel(depth)

def test_fit_point_cloud():
    """Test fit_point_cloud: fitting a set of points to a point cloud"""
    # evenly spaced target points on a sphere
    u = np.linspace(0, np.pi, 150)
    v = np.linspace(0, np.pi, 150)

    x = np.outer(np.cos(u), np.sin(v)).reshape((-1, 1))
    y = np.outer(np.sin(u), np.sin(v)).reshape((-1, 1))
    z = np.outer(np.ones(np.size(u)), np.cos(v)).reshape((-1, 1)) * 3

    tgt_pts = np.hstack((x, y, z))
    tgt_pts = _decimate_points(tgt_pts, .05)

    # pick some points to fit
    some_tgt_pts = tgt_pts[::362]

    # rotation only
    trans = rotation(1.5, .3, -0.4)
    src_pts = apply_trans(trans, some_tgt_pts)
    trans_est = fit_point_cloud(src_pts, tgt_pts, rotate=True, translate=False,
                                scale=0, out='trans')
    est_pts = apply_trans(trans_est, src_pts)
    err = _point_cloud_error(est_pts, tgt_pts)
    assert_array_less(err, .1, "fit_point_cloud with rotation.")

    # rotation and translation
    trans = np.dot(rotation(0.5, .3, -0.4), translation(.3, .2, -.2))
    src_pts = apply_trans(trans, some_tgt_pts)
    trans_est = fit_point_cloud(src_pts, tgt_pts, rotate=True, translate=True,
                                scale=0, out='trans')
    est_pts = apply_trans(trans_est, src_pts)
    err = _point_cloud_error(est_pts, tgt_pts)
    assert_array_less(err, .1, "fit_point_cloud with rotation and "
                      "translation.")

    # rotation and 1 scale parameter
    trans = np.dot(rotation(0.5, .3, -0.4), scaling(1.5, 1.5, 1.5))
    src_pts = apply_trans(trans, some_tgt_pts)
    trans_est = fit_point_cloud(src_pts, tgt_pts, rotate=True, translate=False,
                                scale=1, out='trans')
    est_pts = apply_trans(trans_est, src_pts)
    err = _point_cloud_error(est_pts, tgt_pts)
    assert_array_less(err, .1, "fit_point_cloud with rotation and 1 scaling "
                      "parameter.")

    # rotation and 3 scale parameter
    trans = np.dot(rotation(0.5, .3, -0.4), scaling(1.5, 1.7, 1.1))
    src_pts = apply_trans(trans, some_tgt_pts)
    trans_est = fit_point_cloud(src_pts, tgt_pts, rotate=True, translate=False,
                                scale=3, out='trans')
    est_pts = apply_trans(trans_est, src_pts)
    err = _point_cloud_error(est_pts, tgt_pts)
    assert_array_less(err, .1, "fit_point_cloud with rotation and 3 scaling "
                      "parameters.")

def test_hsp_elp():
    """Test KIT usage of *.elp and *.hsp files against *.txt files."""
    raw_txt = read_raw_kit(sqd_path, mrk_path, elp_txt_path, hsp_txt_path)
    raw_elp = read_raw_kit(sqd_path, mrk_path, elp_path, hsp_path)

    # head points
    pts_txt = np.array([dig_point['r'] for dig_point in raw_txt.info['dig']])
    pts_elp = np.array([dig_point['r'] for dig_point in raw_elp.info['dig']])
    assert_array_almost_equal(pts_elp, pts_txt, decimal=5)

    # transforms
    trans_txt = raw_txt.info['dev_head_t']['trans']
    trans_elp = raw_elp.info['dev_head_t']['trans']
    assert_array_almost_equal(trans_elp, trans_txt, decimal=5)

    # head points in device space
    pts_txt_in_dev = apply_trans(linalg.inv(trans_txt), pts_txt)
    pts_elp_in_dev = apply_trans(linalg.inv(trans_elp), pts_elp)
    assert_array_almost_equal(pts_elp_in_dev, pts_txt_in_dev, decimal=5)

def test_set_dig_montage():
    """Test applying DigMontage to inst."""
    # Extensive testing of applying `dig` to info is done in test_meas_info
    # with `test_make_dig_points`.
    names = ['nasion', 'lpa', 'rpa', '1', '2', '3', '4', '5']
    hsp_points = _read_dig_points(hsp)
    elp_points = _read_dig_points(elp)
    nasion, lpa, rpa = elp_points[:3]
    nm_trans = get_ras_to_neuromag_trans(nasion, lpa, rpa)
    elp_points = apply_trans(nm_trans, elp_points)
    nasion, lpa, rpa = elp_points[:3]
    hsp_points = apply_trans(nm_trans, hsp_points)

    montage = read_dig_montage(hsp, hpi, elp, names, transform=True,
                               dev_head_t=True)
    temp_dir = _TempDir()
    fname_temp = op.join(temp_dir, 'test.fif')
    montage.save(fname_temp)
    montage_read = read_dig_montage(fif=fname_temp)
    for use_mon in (montage, montage_read):
        info = create_info(['Test Ch'], 1e3, ['eeg'])
        with pytest.warns(None):  # warns on one run about not all positions
            _set_montage(info, use_mon)
        hs = np.array([p['r'] for i, p in enumerate(info['dig'])
                       if p['kind'] == FIFF.FIFFV_POINT_EXTRA])
        nasion_dig = np.array([p['r'] for p in info['dig']
                               if all([p['ident'] == FIFF.FIFFV_POINT_NASION,
                                       p['kind'] == FIFF.FIFFV_POINT_CARDINAL])
                               ])
        lpa_dig = np.array([p['r'] for p in info['dig']
                            if all([p['ident'] == FIFF.FIFFV_POINT_LPA,
                                    p['kind'] == FIFF.FIFFV_POINT_CARDINAL])])
        rpa_dig = np.array([p['r'] for p in info['dig']
                            if all([p['ident'] == FIFF.FIFFV_POINT_RPA,
                                    p['kind'] == FIFF.FIFFV_POINT_CARDINAL])])
        hpi_dig = np.array([p['r'] for p in info['dig']
                            if p['kind'] == FIFF.FIFFV_POINT_HPI])
        assert_allclose(hs, hsp_points, atol=1e-7)
        assert_allclose(nasion_dig.ravel(), nasion, atol=1e-7)
        assert_allclose(lpa_dig.ravel(), lpa, atol=1e-7)
        assert_allclose(rpa_dig.ravel(), rpa, atol=1e-7)
        assert_allclose(hpi_dig, elp_points[3:], atol=1e-7)

def test_get_ras_to_neuromag_trans():
    """Test the coordinate transformation from ras to neuromag"""
    # create model points in neuromag-like space
    anterior = [0, 1, 0]
    left = [-1, 0, 0]
    right = [0.8, 0, 0]
    up = [0, 0, 1]
    rand_pts = np.random.uniform(-1, 1, (3, 3))
    pts = np.vstack((anterior, left, right, up, rand_pts))

    # change coord system
    rx, ry, rz, tx, ty, tz = np.random.uniform(-2 * np.pi, 2 * np.pi, 6)
    trans = np.dot(translation(tx, ty, tz), rotation(rx, ry, rz))
    pts_changed = apply_trans(trans, pts)

    # transform back into original space
    nas, lpa, rpa = pts_changed[:3]
    hsp_trans = get_ras_to_neuromag_trans(nas, lpa, rpa)
    pts_restored = apply_trans(hsp_trans, pts_changed)

    err = "Neuromag transformation failed"
    assert_array_almost_equal(pts_restored, pts, 6, err)

def test_set_dig_montage():
    """Test applying DigMontage to inst

    Extensive testing of applying `dig` to info is done in test_meas_info
    with `test_make_dig_points`.
    """
    names = ['nasion', 'lpa', 'rpa', '1', '2', '3', '4', '5']
    hsp_points = _read_dig_points(hsp)
    elp_points = _read_dig_points(elp)
    hpi_points = read_mrk(hpi)
    p0, p1, p2 = elp_points[:3]
    nm_trans = get_ras_to_neuromag_trans(p0, p1, p2)
    elp_points = apply_trans(nm_trans, elp_points)
    nasion_point, lpa_point, rpa_point = elp_points[:3]
    hsp_points = apply_trans(nm_trans, hsp_points)

    montage = read_dig_montage(hsp, hpi, elp, names, unit='m', transform=True)
    info = create_info(['Test Ch'], 1e3, ['eeg'])
    _set_montage(info, montage)
    hs = np.array([p['r'] for i, p in enumerate(info['dig'])
                   if p['kind'] == FIFF.FIFFV_POINT_EXTRA])
    nasion_dig = np.array([p['r'] for p in info['dig']
                           if all([p['ident'] == FIFF.FIFFV_POINT_NASION,
                                   p['kind'] == FIFF.FIFFV_POINT_CARDINAL])])
    lpa_dig = np.array([p['r'] for p in info['dig']
                        if all([p['ident'] == FIFF.FIFFV_POINT_LPA,
                                p['kind'] == FIFF.FIFFV_POINT_CARDINAL])])
    rpa_dig = np.array([p['r'] for p in info['dig']
                        if all([p['ident'] == FIFF.FIFFV_POINT_RPA,
                                p['kind'] == FIFF.FIFFV_POINT_CARDINAL])])
    hpi_dig = np.array([p['r'] for p in info['dig']
                        if p['kind'] == FIFF.FIFFV_POINT_HPI])
    assert_array_equal(hs, hsp_points)
    assert_array_equal(nasion_dig.ravel(), nasion_point)
    assert_array_equal(lpa_dig.ravel(), lpa_point)
    assert_array_equal(rpa_dig.ravel(), rpa_point)
    assert_array_equal(hpi_dig, hpi_points)
    assert_array_equal(montage.dev_head_t, info['dev_head_t']['trans'])

def test_head_to_mni():
    """Test conversion of aseg vertices to MNI coordinates."""
    # obtained using freeview
    coords = np.array([[22.52, 11.24, 17.72], [22.52, 5.46, 21.58],
                       [16.10, 5.46, 22.23], [21.24, 8.36, 22.23]])

    xfm = _read_talxfm('sample', subjects_dir)
    coords_MNI = apply_trans(xfm['trans'], coords)

    trans = read_trans(trans_fname)  # head->MRI (surface RAS)
    mri_head_t = invert_transform(trans)  # MRI (surface RAS)->head matrix

    # obtained from sample_audvis-meg-oct-6-mixed-fwd.fif
    coo_right_amygdala = np.array([[0.01745682,  0.02665809,  0.03281873],
                                   [0.01014125,  0.02496262,  0.04233755],
                                   [0.01713642,  0.02505193,  0.04258181],
                                   [0.01720631,  0.03073877,  0.03850075]])
    coords_MNI_2 = head_to_mni(coo_right_amygdala, 'sample', mri_head_t,
                               subjects_dir)
    # less than 1mm error
    assert_allclose(coords_MNI, coords_MNI_2, atol=10.0)

def test_coregister_fiducials():
    """Test coreg.coregister_fiducials()."""
    # prepare head and MRI fiducials
    trans = Transform('head', 'mri',
                      rotation(.4, .1, 0).dot(translation(.1, -.1, .1)))
    coords_orig = np.array([[-0.08061612, -0.02908875, -0.04131077],
                            [0.00146763, 0.08506715, -0.03483611],
                            [0.08436285, -0.02850276, -0.04127743]])
    coords_trans = apply_trans(trans, coords_orig)

    def make_dig(coords, cf):
        return ({'coord_frame': cf, 'ident': 1, 'kind': 1, 'r': coords[0]},
                {'coord_frame': cf, 'ident': 2, 'kind': 1, 'r': coords[1]},
                {'coord_frame': cf, 'ident': 3, 'kind': 1, 'r': coords[2]})

    mri_fiducials = make_dig(coords_trans, FIFF.FIFFV_COORD_MRI)
    info = {'dig': make_dig(coords_orig, FIFF.FIFFV_COORD_HEAD)}

    # test coregister_fiducials()
    trans_est = coregister_fiducials(info, mri_fiducials)
    assert trans_est.from_str == trans.from_str
    assert trans_est.to_str == trans.to_str
    assert_array_almost_equal(trans_est['trans'], trans['trans'])

def test_read_ctf():
    """Test CTF reader."""
    temp_dir = _TempDir()
    out_fname = op.join(temp_dir, 'test_py_raw.fif')

    # Create a dummy .eeg file so we can test our reading/application of it
    os.mkdir(op.join(temp_dir, 'randpos'))
    ctf_eeg_fname = op.join(temp_dir, 'randpos', ctf_fname_catch)
    shutil.copytree(op.join(ctf_dir, ctf_fname_catch), ctf_eeg_fname)
    with pytest.warns(RuntimeWarning, match='RMSP .* changed to a MISC ch'):
        raw = _test_raw_reader(read_raw_ctf, directory=ctf_eeg_fname)
    picks = pick_types(raw.info, meg=False, eeg=True)
    pos = np.random.RandomState(42).randn(len(picks), 3)
    fake_eeg_fname = op.join(ctf_eeg_fname, 'catch-alp-good-f.eeg')
    # Create a bad file
    with open(fake_eeg_fname, 'wb') as fid:
        fid.write('foo\n'.encode('ascii'))
    pytest.raises(RuntimeError, read_raw_ctf, ctf_eeg_fname)
    # Create a good file
    with open(fake_eeg_fname, 'wb') as fid:
        for ii, ch_num in enumerate(picks):
            args = (str(ch_num + 1), raw.ch_names[ch_num],) + tuple(
                '%0.5f' % x for x in 100 * pos[ii])  # convert to cm
            fid.write(('\t'.join(args) + '\n').encode('ascii'))
    pos_read_old = np.array([raw.info['chs'][p]['loc'][:3] for p in picks])
    with pytest.warns(RuntimeWarning, match='RMSP .* changed to a MISC ch'):
        raw = read_raw_ctf(ctf_eeg_fname)  # read modified data
    pos_read = np.array([raw.info['chs'][p]['loc'][:3] for p in picks])
    assert_allclose(apply_trans(raw.info['ctf_head_t'], pos), pos_read,
                    rtol=1e-5, atol=1e-5)
    assert (pos_read == pos_read_old).mean() < 0.1
    shutil.copy(op.join(ctf_dir, 'catch-alp-good-f.ds_randpos_raw.fif'),
                op.join(temp_dir, 'randpos', 'catch-alp-good-f.ds_raw.fif'))

    # Create a version with no hc, starting out *with* EEG pos (error)
    os.mkdir(op.join(temp_dir, 'nohc'))
    ctf_no_hc_fname = op.join(temp_dir, 'no_hc', ctf_fname_catch)
    shutil.copytree(ctf_eeg_fname, ctf_no_hc_fname)
    remove_base = op.join(ctf_no_hc_fname, op.basename(ctf_fname_catch[:-3]))
    os.remove(remove_base + '.hc')
    with pytest.warns(RuntimeWarning, match='MISC channel'):
        pytest.raises(RuntimeError, read_raw_ctf, ctf_no_hc_fname)
    os.remove(remove_base + '.eeg')
    shutil.copy(op.join(ctf_dir, 'catch-alp-good-f.ds_nohc_raw.fif'),
                op.join(temp_dir, 'no_hc', 'catch-alp-good-f.ds_raw.fif'))

    # All our files
    use_fnames = [op.join(ctf_dir, c) for c in ctf_fnames]
    for fname in use_fnames:
        raw_c = read_raw_fif(fname + '_raw.fif', preload=True)
        with pytest.warns(None):  # sometimes matches "MISC channel"
            raw = read_raw_ctf(fname)

        # check info match
        assert_array_equal(raw.ch_names, raw_c.ch_names)
        assert_allclose(raw.times, raw_c.times)
        assert_allclose(raw._cals, raw_c._cals)
        assert_equal(raw.info['meas_id']['version'],
                     raw_c.info['meas_id']['version'] + 1)
        for t in ('dev_head_t', 'dev_ctf_t', 'ctf_head_t'):
            assert_allclose(raw.info[t]['trans'], raw_c.info[t]['trans'],
                            rtol=1e-4, atol=1e-7)
        for key in ('acq_pars', 'acq_stim', 'bads',
                    'ch_names', 'custom_ref_applied', 'description',
                    'events', 'experimenter', 'highpass', 'line_freq',
                    'lowpass', 'nchan', 'proj_id', 'proj_name',
                    'projs', 'sfreq', 'subject_info'):
            assert_equal(raw.info[key], raw_c.info[key], key)
        if op.basename(fname) not in single_trials:
            # We don't force buffer size to be smaller like MNE-C
            assert raw.buffer_size_sec == raw_c.buffer_size_sec
        assert_equal(len(raw.info['comps']), len(raw_c.info['comps']))
        for c1, c2 in zip(raw.info['comps'], raw_c.info['comps']):
            for key in ('colcals', 'rowcals'):
                assert_allclose(c1[key], c2[key])
            assert_equal(c1['save_calibrated'], c2['save_calibrated'])
            for key in ('row_names', 'col_names', 'nrow', 'ncol'):
                assert_array_equal(c1['data'][key], c2['data'][key])
            assert_allclose(c1['data']['data'], c2['data']['data'], atol=1e-7,
                            rtol=1e-5)
        assert_allclose(raw.info['hpi_results'][0]['coord_trans']['trans'],
                        raw_c.info['hpi_results'][0]['coord_trans']['trans'],
                        rtol=1e-5, atol=1e-7)
        assert_equal(len(raw.info['chs']), len(raw_c.info['chs']))
        for ii, (c1, c2) in enumerate(zip(raw.info['chs'], raw_c.info['chs'])):
            for key in ('kind', 'scanno', 'unit', 'ch_name', 'unit_mul',
                        'range', 'coord_frame', 'coil_type', 'logno'):
                if c1['ch_name'] == 'RMSP' and \
                        'catch-alp-good-f' in fname and \
                        key in ('kind', 'unit', 'coord_frame', 'coil_type',
                                'logno'):
                    continue  # XXX see below...
                assert_equal(c1[key], c2[key], err_msg=key)
            for key in ('cal',):
                assert_allclose(c1[key], c2[key], atol=1e-6, rtol=1e-4,
                                err_msg='raw.info["chs"][%d][%s]' % (ii, key))
            # XXX 2016/02/24: fixed bug with normal computation that used
            # to exist, once mne-C tools are updated we should update our FIF
            # conversion files, then the slices can go away (and the check
            # can be combined with that for "cal")
#.........这里部分代码省略.........

#.........这里部分代码省略.........
                np.linalg.norm, 1,
                montage.pos - centroid)
            assert_array_less(distance_from_centroid, 0.2)
            assert_array_less(0.01, distance_from_centroid)
        assert_array_almost_equal(poss[key], montage.pos, 4, err_msg=key)

    # Test reading in different letter case.
    ch_names = ["F3", "FZ", "F4", "FC3", "FCz", "FC4", "C3", "CZ", "C4", "CP3",
                "CPZ", "CP4", "P3", "PZ", "P4", "O1", "OZ", "O2"]
    montage = read_montage('standard_1020', ch_names=ch_names)
    assert_array_equal(ch_names, montage.ch_names)

    # test transform
    input_strs = ["""
    eeg Fp1 -95.0 -31.0 -3.0
    eeg AF7 -81 -59 -3
    eeg AF3 -87 -41 28
    cardinal 2 -91 0 -42
    cardinal 1 0 -91 -42
    cardinal 3 0 91 -42
    """, """
    Fp1 -95.0 -31.0 -3.0
    AF7 -81 -59 -3
    AF3 -87 -41 28
    nasion -91 0 -42
    lpa 0 -91 -42
    rpa 0 91 -42
    """]

    all_fiducials = [['2', '1', '3'], ['nasion', 'lpa', 'rpa']]

    kinds = ['test_fid.hpts',  'test_fid.sfp']

    for kind, fiducials, input_str in zip(kinds, all_fiducials, input_strs):
        fname = op.join(tempdir, kind)
        with open(fname, 'w') as fid:
            fid.write(input_str)
        montage = read_montage(op.join(tempdir, kind), transform=True)

        # check coordinate transformation
        pos = np.array([-95.0, -31.0, -3.0])
        nasion = np.array([-91, 0, -42])
        lpa = np.array([0, -91, -42])
        rpa = np.array([0, 91, -42])
        fids = np.vstack((nasion, lpa, rpa))
        trans = get_ras_to_neuromag_trans(fids[0], fids[1], fids[2])
        pos = apply_trans(trans, pos)
        assert_array_equal(montage.pos[0], pos)
        idx = montage.ch_names.index(fiducials[0])
        assert_array_equal(montage.pos[idx, [0, 2]], [0, 0])
        idx = montage.ch_names.index(fiducials[1])
        assert_array_equal(montage.pos[idx, [1, 2]], [0, 0])
        idx = montage.ch_names.index(fiducials[2])
        assert_array_equal(montage.pos[idx, [1, 2]], [0, 0])
        pos = np.array([-95.0, -31.0, -3.0])
        montage_fname = op.join(tempdir, kind)
        montage = read_montage(montage_fname, unit='mm')
        assert_array_equal(montage.pos[0], pos * 1e-3)

        # test with last
        info = create_info(montage.ch_names, 1e3,
                           ['eeg'] * len(montage.ch_names))
        _set_montage(info, montage)
        pos2 = np.array([c['loc'][:3] for c in info['chs']])
        assert_array_equal(pos2, montage.pos)
        assert_equal(montage.ch_names, info['ch_names'])

        info = create_info(
            montage.ch_names, 1e3, ['eeg'] * len(montage.ch_names))

        evoked = EvokedArray(
            data=np.zeros((len(montage.ch_names), 1)), info=info, tmin=0)
        evoked.set_montage(montage)
        pos3 = np.array([c['loc'][:3] for c in evoked.info['chs']])
        assert_array_equal(pos3, montage.pos)
        assert_equal(montage.ch_names, evoked.info['ch_names'])

        # Warning should be raised when some EEG are not specified in montage
        with warnings.catch_warnings(record=True) as w:
            info = create_info(montage.ch_names + ['foo', 'bar'], 1e3,
                               ['eeg'] * (len(montage.ch_names) + 2))
            _set_montage(info, montage)
            assert_true(len(w) == 1)

    # Channel names can be treated case insensitive
    with warnings.catch_warnings(record=True) as w:
        info = create_info(['FP1', 'af7', 'AF3'], 1e3, ['eeg'] * 3)
        _set_montage(info, montage)
        assert_true(len(w) == 0)

    # Unless there is a collision in names
    with warnings.catch_warnings(record=True) as w:
        info = create_info(['FP1', 'Fp1', 'AF3'], 1e3, ['eeg'] * 3)
        _set_montage(info, montage)
        assert_true(len(w) == 1)
    with warnings.catch_warnings(record=True) as w:
        montage.ch_names = ['FP1', 'Fp1', 'AF3']
        info = create_info(['fp1', 'AF3'], 1e3, ['eeg', 'eeg'])
        _set_montage(info, montage)
        assert_true(len(w) == 1)

def plot_visualize_mft_sources(fwdmag, stcdata, tmin, tstep,
                               subject, subjects_dir):
    '''
    Plot the MFT sources at time point of peak.
    '''
    print "##### Attempting to plot:"
    # cf. decoding/plot_decoding_spatio_temporal_source.py
    vertices = [s['vertno'] for s in fwdmag['src']]
    if len(vertices) == 1:
        vertices = [fwdmag['src'][0]['vertno'][fwdmag['src'][0]['rr'][fwdmag['src'][0]['vertno']][:, 0] <= -0.],
                    fwdmag['src'][0]['vertno'][fwdmag['src'][0]['rr'][fwdmag['src'][0]['vertno']][:, 0] > -0.]]

    stc_feat = SourceEstimate(stcdata, vertices=vertices,
                              tmin=-0.2, tstep=tstep, subject=subject)
    for hemi in ['lh', 'rh']:
        brain = stc_feat.plot(surface='white', hemi=hemi, subjects_dir=subjects_dir,
                              transparent=True, clim='auto')
        brain.show_view('lateral')
        # use peak getter to move visualization to the time point of the peak
        tmin = 0.095
        tmax = 0.10
        print "Restricting peak search to [%fs, %fs]" % (tmin, tmax)
        if hemi == 'both':
            vertno_max, time_idx = stc_feat.get_peak(hemi='rh', time_as_index=True,
                                                     tmin=tmin, tmax=tmax)
        else:
            vertno_max, time_idx = stc_feat.get_peak(hemi=hemi, time_as_index=True,
                                                     tmin=tmin, tmax=tmax)
        if hemi == 'lh':
            comax = fwdmag['src'][0]['rr'][vertno_max]
            print "hemi=%s: vertno_max=%d, time_idx=%d fwdmag['src'][0]['rr'][vertno_max] = " %\
                  (hemi, vertno_max, time_idx), comax
        elif len(fwdmag['src']) > 1:
            comax = fwdmag['src'][1]['rr'][vertno_max]
            print "hemi=%s: vertno_max=%d, time_idx=%d fwdmag['src'][1]['rr'][vertno_max] = " %\
                  (hemi, vertno_max, time_idx), comax

        print "hemi=%s: setting time_idx=%d" % (hemi, time_idx)
        brain.set_data_time_index(time_idx)
        # draw marker at maximum peaking vertex
        brain.add_foci(vertno_max, coords_as_verts=True, hemi=hemi, color='blue',
                       scale_factor=0.6)
        offsets = np.append([0], [s['nuse'] for s in fwdmag['src']])
        if hemi == 'lh':
            ifoci = [np.nonzero([stcdata[0:offsets[1],time_idx]>=0.25*np.max(stcdata[:,time_idx])][0])]
            vfoci = fwdmag['src'][0]['vertno'][ifoci[0][0]]
            cfoci = fwdmag['src'][0]['rr'][vfoci]
            print "Coords  of %d sel. vfoci: " % cfoci.shape[0]
            print cfoci
            print "vfoci: "
            print vfoci
            print "brain.geo['lh'].coords[vfoci] : "
            print brain.geo['lh'].coords[vfoci]
        elif len(fwdmag['src']) > 1:
            ifoci = [np.nonzero([stcdata[offsets[1]:,time_idx]>=0.25*np.max(stcdata[:,time_idx])][0])]
            vfoci = fwdmag['src'][1]['vertno'][ifoci[0][0]]
            cfoci = fwdmag['src'][1]['rr'][vfoci]
            print "Coords  of %d sel. vfoci: " % cfoci.shape[0]
            print cfoci
            print "vfoci: "
            print vfoci
            print "brain.geo['rh'].coords[vfoci] : "
            print brain.geo['rh'].coords[vfoci]

        mrfoci = np.zeros(cfoci.shape)
        invmri_head_t = invert_transform(fwdmag['info']['mri_head_t'])
        mrfoci = apply_trans(invmri_head_t['trans'],cfoci, move=True)
        print "mrfoci: "
        print mrfoci

        # Just some blops:
        bloblist = np.zeros((300,3))
        for i in xrange(100):
            bloblist[i,0] = float(i)
            bloblist[i+100,1] = float(i)
            bloblist[i+200,2] = float(i)
        mrblobs = apply_trans(invmri_head_t['trans'], bloblist, move=True)
        brain.save_image('testfig_map_%s.png' % hemi)
        brain.close()