def test_vertex_to_mni_fs_nibabel():
    """Test equivalence of vert_to_mni for nibabel and freesurfer
    """
    n_check = 1000
    for subject in ['sample', 'fsaverage']:
        vertices = np.random.randint(0, 100000, n_check)
        hemis = np.random.randint(0, 1, n_check)
        coords = vertex_to_mni(vertices, hemis, subject, mode='nibabel')
        coords_2 = vertex_to_mni(vertices, hemis, subject, mode='freesurfer')
        # less than 0.1 mm error
        assert_allclose(coords, coords_2, atol=0.1)

def test_vertex_to_mni_fs_nibabel():
    """Test equivalence of vert_to_mni for nibabel and freesurfer."""
    n_check = 1000
    subject = 'sample'
    vertices = rng.randint(0, 100000, n_check)
    hemis = rng.randint(0, 1, n_check)
    coords = vertex_to_mni(vertices, hemis, subject, subjects_dir,
                           'nibabel')
    coords_2 = vertex_to_mni(vertices, hemis, subject, subjects_dir,
                             'freesurfer')
    # less than 0.1 mm error
    assert_allclose(coords, coords_2, atol=0.1)

def plot_perm_ttest_results(events_id, inverse_method='dSPM', plot_type='scatter_plot'):
    print('plot_perm_ttest_results')
    all_data = defaultdict(dict)
    fsave_vertices = [np.arange(10242), np.arange(10242)]
    fs_pts = mne.vertex_to_mni(fsave_vertices, [0, 1], 'fsaverage', LOCAL_SUBJECTS_DIR) # 0 for lh
    for cond_id, cond_name, patient, hc, data in patients_hcs_conds_gen(events_id, True, inverse_method):
        all_data[patient][hc] = data[()]
    print(all_data.keys())
    for patient, pat_data in all_data.iteritems():
        print(patient)
        fol = op.join(LOCAL_ROOT_DIR, 'permutation_ttest_results', patient)
        utils.make_dir(fol)
        if op.isfile(op.join(fol, 'perm_ttest_points.npz')):
            d = np.load(op.join(fol, 'perm_ttest_points.npz'))
            if plot_type == 'scatter_plot':
                points, values = d['points'][()], d['values'][()]
            elif plot_type == 'pysurfer':
                vertices, vertives_values = d['vertices'][()], d['vertives_values'][()]
        else:
            points, values, vertices, vertives_values = calc_points(pat_data, fs_pts)
            np.savez(op.join(fol, 'perm_ttest_points'), points=points, values=values, vertices=vertices, vertives_values=vertives_values)
        max_vals = 8 # int(np.percentile([max(v) for v in values.values()], 70))
        print(max_vals)
        fol = op.join(fol, '{}_figures'.format(plot_type))
        utils.make_dir(fol)
        if plot_type == 'scatter_plot':
            scatter_plot_perm_ttest_results(points, values, fs_pts, max_vals, fol)
        elif plot_type == 'pysurfer':
            pysurfer_plot_perm_ttest_results(vertices, vertives_values, max_vals, fol)

def test_vertex_to_mni():
    """Test conversion of vertices to MNI coordinates."""
    # obtained using "tksurfer (sample) (l/r)h white"
    vertices = [100960, 7620, 150549, 96761]
    coords = np.array([[-60.86, -11.18, -3.19], [-36.46, -93.18, -2.36],
                       [-38.00, 50.08, -10.61], [47.14, 8.01, 46.93]])
    hemis = [0, 0, 0, 1]
    coords_2 = vertex_to_mni(vertices, hemis, 'sample', subjects_dir)
    # less than 1mm error
    assert_allclose(coords, coords_2, atol=1.0)

def test_vertex_to_mni():
    """Test conversion of vertices to MNI coordinates
    """
    # obtained using "tksurfer (sample/fsaverage) (l/r)h white"
    vertices = [100960, 7620, 150549, 96761]
    coords_s = np.array([[-60.86, -11.18, -3.19], [-36.46, -93.18, -2.36],
                         [-38.00, 50.08, -10.61], [47.14, 8.01, 46.93]])
    coords_f = np.array([[-41.28, -40.04, 18.20], [-6.05, 49.74, -18.15],
                         [-61.71, -14.55, 20.52], [21.70, -60.84, 25.02]])
    hemis = [0, 0, 0, 1]
    for coords, subject in zip([coords_s, coords_f], ['sample', 'fsaverage']):
        coords_2 = vertex_to_mni(vertices, hemis, subject, subjects_dir)
        # less than 1mm error
        assert_allclose(coords, coords_2, atol=1.0)

def test_vertex_to_mni():
    """Test conversion of vertices to MNI coordinates
    """
    # these were random vertices pulled from "sample" in mne_analyze
    # but mne_analyze won't load the xfm file! So we must use fsaverage,
    # which is sily because the xfm is the identity matrix..
    #vertices = [109445, 82962, 137444]
    #coords = [[-33.3, 11.5, 50.7], [51.8, -15.4, 30.5], [37.6, 38.4, 57.1]]
    #hemi = [0, 1, 1]
    #coords_2 = vertex_to_mni(vertices, hemis, 'sample')
    vertices = [148611, 157229, 95466]
    coords = [[-55.7, -36.6, -9.6], [-48.5, -35.7, -1.1], [44.0, -34.9, -0.9]]
    hemis = [0, 0, 1]  # , 1]
    coords_2 = np.round(vertex_to_mni(vertices, hemis, 'fsaverage'), 1)
    for vi in range(len(vertices)):
        assert_true(coords[vi] == coords_2[vi].tolist())

def test_scale_mri_xfm():
    """Test scale_mri transforms and MRI scaling."""
    # scale fsaverage
    tempdir = _TempDir()
    os.environ['_MNE_FEW_SURFACES'] = 'true'
    fake_home = testing.data_path()
    # add fsaverage
    create_default_subject(subjects_dir=tempdir, fs_home=fake_home,
                           verbose=True)
    # add sample (with few files)
    sample_dir = op.join(tempdir, 'sample')
    os.mkdir(sample_dir)
    os.mkdir(op.join(sample_dir, 'bem'))
    for dirname in ('mri', 'surf'):
        copytree(op.join(fake_home, 'subjects', 'sample', dirname),
                 op.join(sample_dir, dirname))
    subject_to = 'flachkopf'
    spacing = 'oct2'
    for subject_from in ('fsaverage', 'sample'):
        if subject_from == 'fsaverage':
            scale = 1.  # single dim
        else:
            scale = [0.9, 2, .8]  # separate
        src_from_fname = op.join(tempdir, subject_from, 'bem',
                                 '%s-%s-src.fif' % (subject_from, spacing))
        src_from = mne.setup_source_space(
            subject_from, spacing, subjects_dir=tempdir, add_dist=False)
        write_source_spaces(src_from_fname, src_from)
        print(src_from_fname)
        vertices_from = np.concatenate([s['vertno'] for s in src_from])
        assert len(vertices_from) == 36
        hemis = ([0] * len(src_from[0]['vertno']) +
                 [1] * len(src_from[0]['vertno']))
        mni_from = mne.vertex_to_mni(vertices_from, hemis, subject_from,
                                     subjects_dir=tempdir)
        if subject_from == 'fsaverage':  # identity transform
            source_rr = np.concatenate([s['rr'][s['vertno']]
                                        for s in src_from]) * 1e3
            assert_allclose(mni_from, source_rr)
        if subject_from == 'fsaverage':
            overwrite = skip_fiducials = False
        else:
            with pytest.raises(IOError, match='No fiducials file'):
                scale_mri(subject_from, subject_to,  scale,
                          subjects_dir=tempdir)
            skip_fiducials = True
            with pytest.raises(IOError, match='already exists'):
                scale_mri(subject_from, subject_to,  scale,
                          subjects_dir=tempdir, skip_fiducials=skip_fiducials)
            overwrite = True
        scale_mri(subject_from, subject_to, scale, subjects_dir=tempdir,
                  verbose='debug', overwrite=overwrite,
                  skip_fiducials=skip_fiducials)
        if subject_from == 'fsaverage':
            assert _is_mri_subject(subject_to, tempdir), "Scaling failed"
        src_to_fname = op.join(tempdir, subject_to, 'bem',
                               '%s-%s-src.fif' % (subject_to, spacing))
        assert op.exists(src_to_fname), "Source space was not scaled"
        # Check MRI scaling
        fname_mri = op.join(tempdir, subject_to, 'mri', 'T1.mgz')
        assert op.exists(fname_mri), "MRI was not scaled"
        # Check MNI transform
        src = mne.read_source_spaces(src_to_fname)
        vertices = np.concatenate([s['vertno'] for s in src])
        assert_array_equal(vertices, vertices_from)
        mni = mne.vertex_to_mni(vertices, hemis, subject_to,
                                subjects_dir=tempdir)
        assert_allclose(mni, mni_from, atol=1e-3)  # 0.001 mm
    del os.environ['_MNE_FEW_SURFACES']

# load the stc
stc = mne.read_source_estimate(meg_path + '/sample_audvis-meg')

# load the labels
aud_lh = mne.read_label(meg_path + '/labels/Aud-lh.label')
aud_rh = mne.read_label(meg_path + '/labels/Aud-rh.label')

# extract the time course for different labels from the stc
stc_lh = stc.in_label(aud_lh)
stc_rh = stc.in_label(aud_rh)
stc_bh = stc.in_label(aud_lh + aud_rh)

# calculate center of mass and transform to mni coordinates
vtx, _, t_lh = stc_lh.center_of_mass('sample')
mni_lh = mne.vertex_to_mni(vtx, 0, 'sample')[0]
vtx, _, t_rh = stc_rh.center_of_mass('sample')
mni_rh = mne.vertex_to_mni(vtx, 1, 'sample')[0]

# plot the activation
plt.figure()
plt.axes([.1, .275, .85, .625])
hl = plt.plot(stc.times, stc_lh.data.mean(0), 'b')[0]
hr = plt.plot(stc.times, stc_rh.data.mean(0), 'g')[0]
hb = plt.plot(stc.times, stc_bh.data.mean(0), 'r')[0]
plt.xlabel('Time (s)')
plt.ylabel('Source amplitude (dSPM)')
plt.xlim(stc.times[0], stc.times[-1])

# add a legend including center-of-mass mni coordinates to the plot
labels = ['LH: center of mass = %s' % mni_lh.round(2),

def _cluster_rois(sel_path, label_list, stc, src, min_dist, mni_subject='fsaverage'):
    """
    subfunctions of merge_ROIs
    ----------
    mer_path: str
        The directory for storing merged ROIs.
    label_list: list
        Labels to be merged
    """
    class_list = []
    label_list = _sortlist(label_list, stc, src)
    class_list.append(label_list[0])
    for test_fn in label_list[1:]:
        test_label = mne.read_label(test_fn)
        i = 0
        belong = False
        while (i < len(class_list)) and (belong is False):
            class_label = mne.read_label(class_list[i])
            if test_label.hemi != class_label.hemi:
                i = i + 1
                continue
            else:
                # Get the centroids
                class_stc = stc.in_label(class_label)
                test_stc = stc.in_label(test_label)
                class_pca = stc.extract_label_time_course(class_label, src, mode='pca_flip')
                test_pca = stc.extract_label_time_course(test_label, src, mode='pca_flip')
                class_pca = np.squeeze(class_pca)
                test_pca = np.squeeze(test_pca)
                class_pow = np.sum(class_pca ** 2)
                test_pow = np.sum(test_pca ** 2)
                
                if class_label.hemi == 'lh':
                    h = 0
                elif class_label.hemi == 'rh':
                    h = 1
                class_seed, _, _ = class_stc.center_of_mass(mni_subject, hemi=h)
                test_seed, _, _ = test_stc.center_of_mass(mni_subject, hemi=h)
                class_mni = mne.vertex_to_mni(class_seed, h, mni_subject)[0]
                test_mni = mne.vertex_to_mni(test_seed, h, mni_subject)[0]
                
                exch = False
                if np.max(class_pow) < np.max(test_pow):
                    exch = True
                
                # Compute the centroids distance   
                if np.linalg.norm(class_mni - test_mni) < min_dist:
                    if exch == True:
                        os.remove(class_list[i])
                        class_list[i] = test_fn
                        class_list = _sortlist(class_list, stc, src)
                    elif exch == False:
                        os.remove(test_fn)
                    belong = True
                i = i + 1
                
        if belong is False:
            class_list.append(test_fn)
            class_list = _sortlist(class_list, stc, src)
                
    return len(class_list)

def _cluster_rois(sel_path, label_list, stc, src, min_dist, weight, mni_subject='fsaverage'):
    """
    subfunctions of merge_ROIs
    ----------
    mer_path: str
        The directory for storing merged ROIs.
    label_list: list
        Labels to be merged
    """
    class_list = []
    class_list.append(label_list[0])
    for test_fn in label_list[1:]:
        test_label = mne.read_label(test_fn)
        i = 0
        belong = False
        while (i < len(class_list)) and (belong is False):
            class_label = mne.read_label(class_list[i])
            if test_label.hemi != class_label.hemi:
                i = i + 1
                continue
            else:
                # Get the centroids
                class_stc = stc.in_label(class_label)
                test_stc = stc.in_label(test_label)
                class_pow = np.sum(class_stc.data ** 2, axis=1)
                test_pow = np.sum(test_stc.data ** 2, axis=1)
                if class_label.hemi == 'lh':
                    h = 0
                elif class_label.hemi == 'rh':
                    h = 1
                class_seed = class_stc.vertices[h][np.argmax(class_pow)]
                test_seed = test_stc.vertices[h][np.argmax(test_pow)]
                class_vtx = np.searchsorted(stc.vertices[h], class_seed)
                test_vtx = np.searchsorted(stc.vertices[h], test_seed) 
                class_mni = mne.vertex_to_mni(class_vtx, h, mni_subject)[0]
                test_mni = mne.vertex_to_mni(test_vtx, h, mni_subject)[0]
                
                # Get the representative STCs for class label and test label
                class_ts = stc.data[class_vtx, :]
                test_ts = stc.data[test_vtx, :]
                
                # Mark the more apparent ROI
                max_ts = class_ts
                exch = False
                if np.max(class_pow) < np.max(test_pow):
                    max_ts = test_ts
                    exch = True
                
                # Compute the centroids distance   
                if np.linalg.norm(class_mni - test_mni) < min_dist:
                    if exch == True:
                        os.remove(class_list[i])
                        class_list[i] = test_fn
                    elif exch == False:
                        os.remove(test_fn)
                    belong = True
                
                # Compute the similarity
                else:
                    thre = max_ts.std() * weight
                    diff =  np.abs(np.linalg.norm(class_ts) - np.linalg.norm(test_ts))
                    if diff < thre:
                        if exch == True:
                            os.remove(class_list[i])
                            class_list[i] = test_fn
                        elif exch == False:
                            os.remove(test_fn)
                        belong = True
                i = i + 1
                
        if belong is False:
            class_list.append(test_fn)
                
    return len(class_list)

persistent = [line.rstrip() for line in pfid]
remitted = [line.rstrip() for line in rfid]
sx = np.recfromcsv(sx_fname)
hi = [rec[2] for s in subjs for rec in sx if rec[0]==s]
inatt = [rec[1] for s in subjs for rec in sx if rec[0]==s]

# find closest source to seed
if seed[0] < 0:
    hemis = 0  # LH
else:
    hemis = 1  # RH
# important to do this by source, because not every vertex has a source
fname = data_dir + 'morphed-lcmv-%dto%d-'%(bands[0][0],bands[0][1]) + subjs[0]
stc = mne.read_source_estimate(fname)
# this is what we get when we read in the fsaverage subject
coord = mne.vertex_to_mni(vertices=stc.vertno[hemis],hemis=hemis,subject='fsaverage')
dist = np.sqrt((coord[:,0] - seed[0])**2 + (coord[:,1] - seed[1])**2 + (coord[:,2] - seed[2])**2)
seed_src = np.argmin(dist) + hemis*len(stc.lh_vertno)
print 'Distance to seed: %.2fmm'%np.min(dist)

# for each band, compute subject-based correlation map
for l_freq, h_freq in bands:
    nv_corrs = []
    adhd_corrs = []
    per_corrs = []
    rem_corrs = []
    print 'Band %d to %d Hz'%(l_freq, h_freq)
    cnt=0
    for s in subjs:
        print cnt+1, '/', len(subjs)
        fname = data_dir + 'morphed-lcmv-%dto%d-'%(l_freq,h_freq) + s

def _cluster_rois(sel_path, label_list, stc, src, min_dist, weight, mni_subject='fsaverage'):
    """
    subfunctions of merge_ROIs
    ----------
    mer_path: str
        The directory for storing merged ROIs.
    label_list: list
        Labels to be merged
    """
    class_list = []
    class_list.append(label_list[0])
    for test_fn in label_list[1:]:
        test_label = mne.read_label(test_fn)
        i = 0
        belong = False
        while (i < len(class_list)) and (belong is False):
            class_label = mne.read_label(class_list[i])
            if test_label.hemi != class_label.hemi:
                i = i + 1
                continue
            else:
                class_pca = stc.extract_label_time_course(class_label, src, mode='pca_flip')
                test_pca = stc.extract_label_time_course(test_label, src, mode='pca_flip')
                class_pow = np.sum(class_pca ** 2)
                test_pow = np.sum(test_pca ** 2)
                max_pca = class_pca
                exch = False
                if class_pow < test_pow:
                    max_pca = test_pca
                    exch = True
                    
                nearby = False
                class_stc = stc.in_label(class_label)
                test_stc = stc.in_label(test_label)
                if class_label.hemi == 'lh':
                    class_vtx, _ = class_stc.get_peak(hemi='lh')
                    test_vtx, _ = test_stc.get_peak(hemi='lh') 
                    class_mni = mne.vertex_to_mni(class_vtx, 0, mni_subject)[0]
                    test_mni = mne.vertex_to_mni(test_vtx, 0, mni_subject)[0]
                elif class_label.hemi == 'rh':
                    class_vtx, _ = class_stc.get_peak(hemi='rh')
                    test_vtx, _ = test_stc.get_peak(hemi='rh') 
                    class_mni = mne.vertex_to_mni(class_vtx, 1, mni_subject)[0]
                    test_mni = mne.vertex_to_mni(test_vtx, 1, mni_subject)[0]
                if np.linalg.norm(class_mni - test_mni) < min_dist:
                    if exch == True:
                        os.remove(class_list[i])
                        class_list[i] = test_fn
                    elif exch == False:
                        os.remove(test_fn)
                    nearby = True
                    belong = True
                    
                if nearby == False:
                    thre = max_pca.std() * weight
                    diff =  np.abs(np.linalg.norm(class_pca) - np.linalg.norm(test_pca))
                    if diff < thre:
                        if exch == True:
                            os.remove(class_list[i])
                            class_list[i] = test_fn
                        elif exch == False:
                            os.remove(test_fn)
                        belong = True
                i = i + 1
        if belong is False:
            class_list.append(test_fn)
                
    return len(class_list)

doMorph = True

subjectDir = os.environ['SUBJECTS_DIR']
modelSubj = 'fsaverage'

# Distance from fMRI-based central point to include in hand motor label
labelRadii = [10]
subject_set = [str(subj) for subj in subj_nums_hcp]

# Load list of labels and pull out pre-central gyrus
labelList = mne.read_labels_from_annot(subject=modelSubj, parc='aparc.a2009s',
                                       hemi='both')
primaryMotor = [l for l in labelList if 'G_precentral-' in l.name]

# Get MNI coords for vertices in L/R hemis
primaryMotorMNI_pos = [mne.vertex_to_mni(primaryMotor[0].vertices, 0, modelSubj),
                       mne.vertex_to_mni(primaryMotor[1].vertices, 1, modelSubj)]

# Find closest point in fsaverage brain space according to
# Witt 2008, Functional neuroimaging correlates of finger-tapping task
# variations coords in Talairach space
# list is [left hemi, right hemi]
hand_knob_pos = [np.atleast_2d(np.array([-38, -26, 50])),
                 np.atleast_2d(np.array([36, -22, 54]))]
SMA = [-4, -8, 52]

dists = []
for orig_label, label_pos, knob_pos in zip(primaryMotor, primaryMotorMNI_pos, hand_knob_pos):
    # Find dist between MNI point and all vertices
    dists.append(np.squeeze(cdist(knob_pos, label_pos, 'euclidean')))
    # Find min dist and index

fid.close()

# prepare the fsaverage subject if we're dealing with AFNI mask
if roi.find('+tlrc') > 0:
    mask_fname = '/lscratch/' + os.environ['SLURM_JOBID'] + '/mask.txt'
    os.system('3dmaskdump -xyz -o ' + mask_fname + ' ' + roi)
    a = np.genfromtxt(mask_fname)
    os.system('rm ' + mask_fname)
    gv = a[a[:, 6] > 0, 3:6]  # separate good voxels in the mask
    # change AFNI results from RAI to LPI
    gv[:, 0] = gv[:, 0] * -1
    gv[:, 1] = gv[:, 1] * -1

    src = mne.setup_source_space(subject='fsaverage_mne', fname=None, spacing='ico5', surface='inflated', n_jobs=2)
    # get left and right coordinates for all the sources
    coord0 = mne.vertex_to_mni(vertices=src[0]['vertno'], hemis=0, subject='fsaverage_mne')
    coord1 = mne.vertex_to_mni(vertices=src[1]['vertno'], hemis=1, subject='fsaverage_mne')
    coord = np.vstack([coord0, coord1])
    # store the index of the sources within min_dist of the mask voxels
    b = []
    for i in range(gv.shape[0]):
        dist = np.sqrt((coord[:, 0] - gv[i, 0]) ** 2 + (coord[:, 1] - gv[i, 1]) ** 2 + (coord[:, 2] - gv[i, 2]) ** 2)
        if min(dist) <= min_dist:
            b.append(np.argmin(dist))
    # create a stc with 1s for the near sources
    d = np.zeros([coord.shape[0], 1])
    d[b] = 1
    stc = mne.SourceEstimate(d, vertices=[src[0]['vertno'], src[1]['vertno']],
                             tmin=0, tstep=1, subject='fsaverage_mne')
    # convert the stc to a label so we can morph it per subject later
    avg_label = mne.stc_to_label(stc, src=src, smooth=True, connected=False)

def _cluster_sel(sel_path, label_list, stc, src, min_dist, weight, mni_subject='fsaverage'):
    """
    subfunctions of sel_ROIs
    ----------
    sel_path: string or list
        The directory for storing selected ROIs.
    label_list: list
        Labels to be selected.
    stc: the object of source estimates.
    src: the object of the common source space
    min_dist: int (mm)
        Least distance between ROIs candidates.
    weight: float
        Euclidean_norms weight related to the larger candidate's standard deviation.
    """
    class_list = []
    class_list.append(label_list[0])
    for test_fn in label_list[1:]:
        test_label = mne.read_label(test_fn)
        i = 0
        belong = False
        while (i < len(class_list)) and (belong is False):
            class_label = mne.read_label(class_list[i])
            if test_label.hemi != class_label.hemi:
                i = i + 1
                continue
            else:
                class_pca = stc.extract_label_time_course(class_label, src, mode='pca_flip')
                test_pca = stc.extract_label_time_course(test_label, src, mode='pca_flip')
                class_pow = np.sum(class_pca ** 2)
                test_pow = np.sum(test_pca ** 2)
                max_pca = class_pca
                exch = False
                if class_pow < test_pow:
                    max_pca = test_pca
                    exch = True
                    
                nearby = False
                class_stc = stc.in_label(class_label)
                test_stc = stc.in_label(test_label)
                if class_label.hemi == 'lh':
                    class_vtx, _ = class_stc.get_peak(hemi='lh')
                    test_vtx, _ = test_stc.get_peak(hemi='lh') 
                    class_mni = mne.vertex_to_mni(class_vtx, 0, mni_subject)[0]
                    test_mni = mne.vertex_to_mni(test_vtx, 0, mni_subject)[0]
                elif class_label.hemi == 'rh':
                    class_vtx, _ = class_stc.get_peak(hemi='rh')
                    test_vtx, _ = test_stc.get_peak(hemi='rh') 
                    class_mni = mne.vertex_to_mni(class_vtx, 1, mni_subject)[0]
                    test_mni = mne.vertex_to_mni(test_vtx, 1, mni_subject)[0]
                if np.linalg.norm(class_mni - test_mni) < min_dist:
                    if exch == True:
                        os.remove(class_list[i])
                        class_list[i] = test_fn
                    elif exch == False:
                        os.remove(test_fn)
                    nearby = True
                    belong = True
                    
                if nearby == False:
                    thre = max_pca.std() * weight
                    diff =  np.abs(np.linalg.norm(class_pca) - np.linalg.norm(test_pca))
                    if diff < thre:
                        if exch == True:
                            os.remove(class_list[i])
                            class_list[i] = test_fn
                        elif exch == False:
                            os.remove(test_fn)
                        belong = True
                i = i + 1
        if belong is False:
            class_list.append(test_fn)
                
    return len(class_list)

subject = 'fsaverage'

# Plot brain in 3D with PySurfer if available
brain = new_stc.plot(subject, hemi='lh', subjects_dir=subjects_dir)
brain.show_view('lateral')

# use peak getter to move vizualization to the time point of the peak
vertno_max, time_idx = new_stc.get_peak(hemi='lh', time_as_index=True)

brain.set_data_time_index(time_idx)

# draw marker at maximum peaking vertex
brain.add_foci(vertno_max, coords_as_verts=True, hemi='lh', color='blue',
               scale_factor=0.6, map_surface='white')

mni_coords = mne.vertex_to_mni(vertno_max, hemis=0, subject=subject,
                               subjects_dir=subjects_dir)
print 'The MNI coords are ', mni_coords

#my_trans = mne.read_trans(?)
#src_pts = apply_trans(trans, some_tgt_pts)

from surfer import utils
utils.coord_to_label(subject, mni_coords[0], label='mycoord',
                     hemi='lh', n_steps=25, map_surface="white")
brain.add_label('mycoord-lh.label', color="darkseagreen", alpha=.8)

# if the new mni_coords are computed
brain.add_foci(mni_coords[0], coords_as_verts=False, hemi='lh', color='red',
               map_surface='white', scale_factor=0.6)