def test_coord_transform_trivial():
    sform = np.eye(4)
    x = np.random.random((10,))
    y = np.random.random((10,))
    z = np.random.random((10,))

    x_, y_, z_ = coord_transform(x, y, z, sform)
    np.testing.assert_array_equal(x, x_)
    np.testing.assert_array_equal(y, y_)
    np.testing.assert_array_equal(z, z_)

    sform[:, -1] = 1
    x_, y_, z_ = coord_transform(x, y, z, sform)
    np.testing.assert_array_equal(x + 1, x_)
    np.testing.assert_array_equal(y + 1, y_)
    np.testing.assert_array_equal(z + 1, z_)

    # Test the output in case of one item array
    x, y, z = x[:1], y[:1], z[:1]
    x_, y_, z_ = coord_transform(x, y, z, sform)
    np.testing.assert_array_equal(x + 1, x_)
    np.testing.assert_array_equal(y + 1, y_)
    np.testing.assert_array_equal(z + 1, z_)

    # Test the output in case of simple items
    x, y, z = x[0], y[0], z[0]
    x_, y_, z_ = coord_transform(x, y, z, sform)
    np.testing.assert_array_equal(x + 1, x_)
    np.testing.assert_array_equal(y + 1, y_)
    np.testing.assert_array_equal(z + 1, z_)

def test_sample_locations():
    # check positions of samples on toy example, with an affine != identity
    # flat horizontal mesh
    mesh = flat_mesh(5, 7)
    affine = np.diagflat([10, 20, 30, 1])
    inv_affine = np.linalg.inv(affine)
    # transform vertices to world space
    vertices = np.asarray(
        resampling.coord_transform(*mesh[0].T, affine=affine)).T
    # compute by hand the true offsets in voxel space
    # (transformed by affine^-1)
    ball_offsets = surface._load_uniform_ball_cloud(10)
    ball_offsets = np.asarray(
        resampling.coord_transform(*ball_offsets.T, affine=inv_affine)).T
    line_offsets = np.zeros((10, 3))
    line_offsets[:, 2] = np.linspace(-1, 1, 10)
    line_offsets = np.asarray(
        resampling.coord_transform(*line_offsets.T, affine=inv_affine)).T
    # check we get the same locations
    for kind, offsets in [('line', line_offsets), ('ball', ball_offsets)]:
        locations = surface._sample_locations(
            [vertices, mesh[1]], affine, 1., kind=kind, n_points=10)
        true_locations = np.asarray([vertex + offsets for vertex in mesh[0]])
        assert_array_equal(locations.shape, true_locations.shape)
        assert_array_almost_equal(true_locations, locations)
    assert_raises(ValueError, surface._sample_locations,
                  mesh, affine, 1., kind='bad_kind')

def get_roi_center(roi_native_path, roi_mni_path):
    """Get ROI center of mass.
    Get back coordinate in img space and in coordinate space.
    Also actual center of mass.
    """
    # computations in native space
    if type(roi_native_path) is str:
        img = nib.load(roi_native_path)
    else:
        img = roi_native_path
    data = img.get_data()
    data = as_ndarray(data)
    my_map = data.copy()
    center_coords = ndimage.center_of_mass(np.abs(my_map))

    x_map, y_map, z_map = center_coords[:3]
    native_coords = np.asarray(coord_transform(x_map, y_map, z_map,
                                               img.get_affine())).tolist()
    voxel = [round(x) for x in center_coords]
    # computations in mni space
    if type(roi_mni_path) is str:
        img = nib.load(roi_mni_path)
    else:
        img = roi_mni_path
    data = img.get_data()
    data = as_ndarray(data)
    my_map = data.copy()
    mni_center_coords = ndimage.center_of_mass(np.abs(my_map))
    x_map, y_map, z_map = mni_center_coords[:3]
    mni_coords = np.asarray(coord_transform(x_map, y_map, z_map,
                                            img.get_affine())).tolist()
    # returns voxel and true center mass coords
    # returns also native and mni space coords
    return (voxel[:3], center_coords[:3], [round(x) for x in native_coords],
            [round(x) for x in mni_coords])

def _get_affinity(seeds, coords, radius, allow_overlap, affine, mask_img=None):
    
    seeds = list(seeds)

    # Compute world coordinates of all in-mask voxels.           
    mask_coords = list(zip(*coords.T))
    # For each seed, get coordinates of nearest voxel
    nearests = []
    for sx, sy, sz in seeds:
        nearest = np.round(coord_transform(sx, sy, sz, np.linalg.inv(affine)))
        nearest = nearest.astype(int)
        nearest = (nearest[0], nearest[1], nearest[2])
        try:
            nearests.append(mask_coords.index(nearest))
        except ValueError:
            nearests.append(None)

    mask_coords = np.asarray(list(zip(*mask_coords)))
    mask_coords = coord_transform(mask_coords[0], mask_coords[1],
                                  mask_coords[2], affine)
    mask_coords = np.asarray(mask_coords).T

    if (radius is not None and
            LooseVersion(sklearn.__version__) < LooseVersion('0.16')):
        # Fix for scikit learn versions below 0.16. See
        # https://github.com/scikit-learn/scikit-learn/issues/4072
        radius += 1e-6

    clf = neighbors.NearestNeighbors(radius=radius)
    A = clf.fit(mask_coords).radius_neighbors_graph(seeds)
    A = A.tolil()
    for i, nearest in enumerate(nearests):
        if nearest is None:
            continue
        A[i, nearest] = True

    # Include the voxel containing the seed itself if not masked
    mask_coords = mask_coords.astype(int).tolist()
    for i, seed in enumerate(seeds):
        try:
            A[i, mask_coords.index(seed)] = True
        except ValueError:
            # seed is not in the mask
            pass

    if not allow_overlap:
        if np.any(A.sum(axis=0) >= 2):
            raise ValueError('Overlap detected between spheres')

    return A

def test_coord_transform_trivial():
    sform = np.eye(4)
    x = np.random.random((10,))
    y = np.random.random((10,))
    z = np.random.random((10,))

    x_, y_, z_ = coord_transform(x, y, z, sform)
    np.testing.assert_array_equal(x, x_)
    np.testing.assert_array_equal(y, y_)
    np.testing.assert_array_equal(z, z_)

    sform[:, -1] = 1
    x_, y_, z_ = coord_transform(x, y, z, sform)
    np.testing.assert_array_equal(x + 1, x_)
    np.testing.assert_array_equal(y + 1, y_)
    np.testing.assert_array_equal(z + 1, z_)

def demo_plot_roi(**kwargs):
    """ Demo plotting an ROI
    """
    mni_affine = MNI152TEMPLATE.get_affine()
    data = np.zeros((91, 109, 91))
    # Color a asymetric rectangle around Broca area:
    x, y, z = -52, 10, 22
    x_map, y_map, z_map = coord_transform(x, y, z, np.linalg.inv(mni_affine))
    data[int(x_map) - 5 : int(x_map) + 5, int(y_map) - 3 : int(y_map) + 3, int(z_map) - 10 : int(z_map) + 10] = 1
    img = nibabel.Nifti1Image(data, mni_affine)
    return plot_roi(img, title="Broca's area", **kwargs)

def get_seeds(ds, radius):
    
    if check_proximity(ds, radius):
        return load_proximity(ds, radius)
    
    
    # Get the seeds
    process_mask_coords = ds.fa.voxel_indices.T
    process_mask_coords = coord_transform(
        process_mask_coords[0], process_mask_coords[1],
        process_mask_coords[2], ds.a.imgaffine)
    process_mask_coords = np.asarray(process_mask_coords).T
    
    seeds = process_mask_coords
    coords = ds.fa.voxel_indices
    logger.info("Building proximity matrix...")
    A = _get_affinity(seeds, coords, radius, allow_overlap=True, affine=ds.a.imgaffine)
    
    save_proximity(ds, radius, A)
    
    return A

def test_outlier_cut_coords():
    """ Test to plot a subset of a large set of cuts found for a small area."""
    bg_img = load_mni152_template()

    data = np.zeros((79, 95, 79))
    affine = np.array([[  -2.,    0.,    0.,   78.],
                       [   0.,    2.,    0., -112.],
                       [   0.,    0.,    2.,  -70.],
                       [   0.,    0.,    0.,    1.]])

    # Color a cube around a corner area:
    x, y, z = 20, 22, 60
    x_map, y_map, z_map = coord_transform(x, y, z,
                                          np.linalg.inv(affine))

    data[int(x_map) - 1:int(x_map) + 1,
         int(y_map) - 1:int(y_map) + 1,
         int(z_map) - 1:int(z_map) + 1] = 1
    img = nibabel.Nifti1Image(data, affine)
    cuts = find_cut_slices(img, n_cuts=20, direction='z')

    p = plot_stat_map(img, display_mode='z', cut_coords=cuts[-4:],
                      bg_img=bg_img)

neg_log_pvals_permuted_ols, _, _ = permuted_ols(
    tested_var, fmri_masked,
    model_intercept=True,
    n_perm=5000,  # 5,000 for the sake of time. Idealy, this should be 10,000
    n_jobs=1)  # can be changed to use more CPUs
neg_log_pvals_permuted_ols_unmasked = nifti_masker.inverse_transform(
    np.ravel(neg_log_pvals_permuted_ols))

### Visualization #############################################################
from nilearn.plotting import plot_stat_map

# Various plotting parameters
z_slice = 12  # plotted slice
from nilearn.image.resampling import coord_transform
affine = neg_log_pvals_anova_unmasked.get_affine()
_, _, k_slice = coord_transform(0, 0, z_slice,
                                linalg.inv(affine))
k_slice = round(k_slice)

threshold = - np.log10(0.1)  # 10% corrected
vmax = min(np.amax(neg_log_pvals_permuted_ols),
           np.amax(neg_log_pvals_anova))

# Plot Anova p-values
fig = plt.figure(figsize=(5, 7), facecolor='k')

display = plot_stat_map(neg_log_pvals_anova_unmasked,
                        threshold=threshold, cmap=plt.cm.autumn,
                        display_mode='z', cut_coords=[z_slice],
                        figure=fig, vmax=vmax, black_bg=True)

neg_log_pvals_anova_data = neg_log_pvals_anova_unmasked.get_data()

def index_to_xy_coord(x, y, z=10):
    '''Transforms data index to coordinates of the background + offset'''
    coords = coord_transform(x, y, z,
                             affine=thresholded_score_map_img.get_affine())
    return np.array(coords)[np.newaxis, :] + np.array([0, 1, 0])

def get_clusters_table(stat_img, stat_threshold, cluster_threshold=None,
                       min_distance=8.):
    """Creates pandas dataframe with img cluster statistics.

    Parameters
    ----------
    stat_img : Niimg-like object,
       Statistical image (presumably in z- or p-scale).

    stat_threshold: `float`
        Cluster forming threshold in same scale as `stat_img` (either a
        p-value or z-scale value).

    cluster_threshold : `int` or `None`, optional
        Cluster size threshold, in voxels.

    min_distance: `float`, optional
        Minimum distance between subpeaks in mm. Default is 8 mm.

    Returns
    -------
    df : `pandas.DataFrame`
        Table with peaks and subpeaks from thresholded `stat_img`. For binary
        clusters (clusters with >1 voxel containing only one value), the table
        reports the center of mass of the cluster, rather than any peaks/subpeaks.
    """
    cols = ['Cluster ID', 'X', 'Y', 'Z', 'Peak Stat', 'Cluster Size (mm3)']
    stat_map = stat_img.get_data()
    conn_mat = np.zeros((3, 3, 3), int)  # 6-connectivity, aka NN1 or "faces"
    conn_mat[1, 1, :] = 1
    conn_mat[1, :, 1] = 1
    conn_mat[:, 1, 1] = 1
    voxel_size = np.prod(stat_img.header.get_zooms())
    
    # Binarize using CDT
    binarized = stat_map > stat_threshold
    binarized = binarized.astype(int)
    
    # If the stat threshold is too high simply return an empty dataframe
    if np.sum(binarized) == 0:
        warnings.warn('Attention: No clusters with stat higher than %f' %
                      stat_threshold)
        return pd.DataFrame(columns=cols)
    
    # Extract connected components above cluster size threshold
    label_map = ndimage.measurements.label(binarized, conn_mat)[0]
    clust_ids = sorted(list(np.unique(label_map)[1:]))
    for c_val in clust_ids:
        if cluster_threshold is not None and np.sum(
                label_map == c_val) < cluster_threshold:
            stat_map[label_map == c_val] = 0
            binarized[label_map == c_val] = 0
    
    # If the cluster threshold is too high simply return an empty dataframe
    # this checks for stats higher than threshold after small clusters
    # were removed from stat_map
    if np.sum(stat_map > stat_threshold) == 0:
        warnings.warn('Attention: No clusters with more than %d voxels' %
                      cluster_threshold)
        return pd.DataFrame(columns=cols)
    
    # Now re-label and create table
    label_map = ndimage.measurements.label(binarized, conn_mat)[0]
    clust_ids = sorted(list(np.unique(label_map)[1:]))
    peak_vals = np.array(
            [np.max(stat_map * (label_map == c)) for c in clust_ids])
    clust_ids = [clust_ids[c] for c in
                 (-peak_vals).argsort()]  # Sort by descending max value
    
    rows = []
    for c_id, c_val in enumerate(clust_ids):
        cluster_mask = label_map == c_val
        masked_data = stat_map * cluster_mask
        
        cluster_size_mm = int(np.sum(cluster_mask) * voxel_size)
        
        # Get peaks, subpeaks and associated statistics
        subpeak_ijk, subpeak_vals = _local_max(masked_data, stat_img.affine,
                                               min_distance=min_distance)
        subpeak_xyz = np.asarray(coord_transform(subpeak_ijk[:, 0],
                                                 subpeak_ijk[:, 1],
                                                 subpeak_ijk[:, 2],
                                                 stat_img.affine)).tolist()
        subpeak_xyz = np.array(subpeak_xyz).T
        
        # Only report peak and, at most, top 3 subpeaks.
        n_subpeaks = np.min((len(subpeak_vals), 4))
        for subpeak in range(n_subpeaks):
            if subpeak == 0:
                row = [c_id + 1, subpeak_xyz[subpeak, 0],
                       subpeak_xyz[subpeak, 1], subpeak_xyz[subpeak, 2],
                       subpeak_vals[subpeak], cluster_size_mm]
            else:
                # Subpeak naming convention is cluster num + letter (1a, 1b, etc.)
                sp_id = '{0}{1}'.format(c_id + 1, ascii_lowercase[subpeak - 1])
                row = [sp_id, subpeak_xyz[subpeak, 0], subpeak_xyz[subpeak, 1],
                       subpeak_xyz[subpeak, 2], subpeak_vals[subpeak], '']
            rows += [row]
    df = pd.DataFrame(columns=cols, data=rows)
    return df

def cluster_stats(stat_img, mask_img, threshold, height_control='fpr',
                  cluster_th=0, nulls=None):
    """
    Return a list of clusters, each cluster being represented by a
    dictionary. Clusters are sorted by descending size order. Within
    each cluster, local maxima are sorted by descending statical value

    Parameters
    ----------
    stat_img: Niimg-like object,
       statsitical image (presumably in z scale)
    mask_img: Niimg-like object,
        mask image
    threshold: float,
        cluster forming threshold (either a p-value or z-scale value)
    height_control: string
        false positive control meaning of cluster forming
        threshold: 'fpr'|'fdr'|'bonferroni'|'none'
    cluster_th: int or float,
        cluster size threshold
    nulls: dictionary,
        statistics of the null distribution

    Notes
    -----
    If there is no cluster, an empty list is returned
    """
    if nulls is None: nulls = {}

    # Masking
    mask_img, stat_img = check_niimg(mask_img), check_niimg(stat_img)
    if not _check_same_fov(mask_img, stat_img):
        raise ValueError('mask_img and stat_img do not have the same fov')
    mask = mask_img.get_data().astype(np.bool)
    affine = mask_img.get_affine()
    stat_map = stat_img.get_data() * mask
    n_voxels = mask.sum()

    # Thresholding
    if height_control == 'fpr':
        z_th = norm.isf(threshold)
    elif height_control == 'fdr':
        z_th = fdr_threshold(stat_map[mask], threshold)
    elif height_control == 'bonferroni':
        z_th = norm.isf(threshold / n_voxels)
    else:  # Brute-force thresholding
        z_th = threshold

    p_th = norm.sf(z_th)
    # General info
    info = {'n_voxels': n_voxels,
            'threshold_z': z_th,
            'threshold_p': p_th,
            'threshold_pcorr': np.minimum(1, p_th * n_voxels)}

    above_th = stat_map > z_th
    above_values = stat_map * above_th
    if (above_th == 0).all():
        return [], info

    # Extract connected components above threshold
    labels, n_labels = label(above_th)

    # Extract the local maxima anove the threshold
    maxima_mask = (above_values ==
                   np.maximum(z_th, maximum_filter(above_values, 3)))
    x, y, z = np.array(np.where(maxima_mask))
    maxima_coords = np.array(coord_transform(x, y, z, affine)).T
    maxima_labels = labels[maxima_mask]
    maxima_values = above_values[maxima_mask]

    # FDR-corrected p-values
    max_fdr_p_values = fdr_p_values(stat_map[mask])[maxima_mask[mask]]

    # Default "nulls"
    if not 'zmax' in nulls:
        nulls['zmax'] = 'bonferroni'
    if not 'smax' in nulls:
        nulls['smax'] = None
    if not 's' in nulls:
        nulls['s'] = None

    # Make list of clusters, each cluster being a dictionary
    clusters = []
    for k in range(n_labels):
        cluster_size = np.sum(labels == k + 1)
        if cluster_size >= cluster_th:

            # get the position of the maxima that belong to that cluster
            in_cluster = maxima_labels == k + 1

            # sort the maxima by decreasing statistical value
            max_vals = maxima_values[in_cluster]
            sorted_ = max_vals.argsort()[::-1]

            # Report significance levels in each cluster
            z_score = max_vals[sorted_]
            p_values = norm.sf(z_score)

            # Voxel-level corrected p-values
#.........这里部分代码省略.........

def coord_transform_z(z, img):
    x, y, z = coord_transform(0, 0, z, img.affine)
    return z

def coord_transform_y(y, img):
    x, y, z = coord_transform(0, y, 0, img.affine)
    return y

def coord_transform_x(x, img):
    x, y, z = coord_transform(x, 0, 0, img.affine)
    return x

def find_region_names_using_cut_coords(coords, atlas_img, labels=None):
    """Given list of MNI space coordinates, get names of the brain regions.

    Names of the brain regions are returned by getting nearest coordinates
    in the given `atlas_img` space iterated over the provided list of
    `coords`. These new image coordinates are then used to grab the label
    number (int) and name assigned to it. Last, these names are returned.

    Parameters
    ----------
    coords : Tuples of coordinates in a list
        MNI coordinates.

    atlas_img : Nifti-like image
        Path to or Nifti-like object. The labels (integers) ordered in
        this image should be sequential. Example: [0, 1, 2, 3, 4] but not
        [0, 5, 6, 7]. Helps in returning correct names without errors.

    labels : str in a list
        Names of the brain regions assigned to each label in atlas_img.
        NOTE: label with index 0 is assumed as background. Example:
            harvard oxford atlas. Hence be removed.

    Returns
    -------
    new_labels : int in a list
        Labels in integers generated according to correspondence with
        given atlas image and provided coordinates.

    names : str in a list
        Names of the brain regions generated according to given inputs.
    """
    if not isinstance(coords, collections.Iterable):
        raise ValueError("coords given must be a list of triplets of "
                         "coordinates in native space [(1, 2, 3)]. "
                         "You provided {0}".format(type(coords)))

    if isinstance(atlas_img, _basestring):
        atlas_img = check_niimg(atlas_img)

    affine = get_affine(atlas_img)
    atlas_data = _safe_get_data(atlas_img, ensure_finite=True)
    check_labels_from_atlas = np.unique(atlas_data)

    if labels is not None:
        names = []
        if not isinstance(labels, collections.Iterable):
            labels = np.asarray(labels)

    if isinstance(labels, collections.Iterable) and \
            isinstance(check_labels_from_atlas, collections.Iterable):
        if len(check_labels_from_atlas) != len(labels):
            warnings.warn("The number of labels provided does not match "
                          "with number of unique labels with atlas image.",
                          stacklevel=2)

    coords = list(coords)
    nearest_coordinates = []

    for sx, sy, sz in coords:
        nearest = np.round(coord_transform(sx, sy, sz, np.linalg.inv(affine)))
        nearest = nearest.astype(int)
        nearest = (nearest[0], nearest[1], nearest[2])
        nearest_coordinates.append(nearest)

    assert(len(nearest_coordinates) == len(coords))

    new_labels = []
    for coord_ in nearest_coordinates:
        # Grab index of current coordinate
        index = atlas_data[coord_]
        new_labels.append(index)
        if labels is not None:
            names.append(labels[index])

    if labels is not None:
        return new_labels, names
    else:
        return new_labels