[MRG] Cortex surface projections

examples/01_plotting/brainpedia_surface.py

@@ -0,0 +1,17 @@
+from matplotlib import pyplot as plt
+
+import nilearn.image
+import nilearn.datasets
+from nilearn.plotting import surf_plotting
+
+brainpedia = nilearn.datasets.fetch_neurovault_ids(image_ids=(32015,))
+image = brainpedia.images[0]
+
+fsaverage = nilearn.datasets.fetch_surf_fsaverage5()
+pial_left = surf_plotting.load_surf_mesh(fsaverage.pial_left)[0]
+
+texture = surf_plotting.niimg_to_surf_data(image, pial_left)
+
+surf_plotting.plot_surf_stat_map(fsaverage.infl_left, texture, cmap='bwr')
+
+plt.show()

examples/01_plotting/show_ball_sampling.py

@@ -0,0 +1,66 @@
+import numpy as np
+import matplotlib.patches
+from matplotlib import pyplot as plt
+from nilearn.plotting import surf_plotting
+
+
+def _piecewise_const_image():
+    w, h = 7, 10
+    image = np.random.uniform(0, 1, size=(w, h))
+    image -= image.min()
+    image /= image.max()
+    big_image = np.empty((10 * w, 10 * h))
+    for i in range(w):
+        for j in range(h):
+            big_image[10 * i:10 * i + 10, 10 * j:10 * j + 10] = image[i, j]
+    return big_image
+
+
+def _random_mesh(image_shape, n_nodes=5):
+    x = np.random.uniform(0, image_shape[0], size=(n_nodes, ))
+    y = np.random.uniform(0, image_shape[1], size=(n_nodes, ))
+    return np.asarray([x, y]).T
+
+
+def show_sampling(ball_radius=10, n_nodes=5, n_points=7, link_points=True):
+    image = _piecewise_const_image()
+    images = [image, image > .5, image < .5, image < .8]
+    mesh = _random_mesh(image.shape, n_nodes=n_nodes)
+    all_values, sample_points, _ = surf_plotting._ball_sampling(
+        images, mesh, np.eye(3), n_points=n_points, ball_radius=ball_radius)
+    fig, axes = plt.subplots(2, 2)
+    axes = axes.ravel()
+    for image, values, ax in zip(images, all_values, axes):
+        ax.imshow(image, cmap='gray', vmin=0, vmax=1)
+        ax.add_patch(
+            matplotlib.patches.Rectangle(
+                (-.5, -.5),
+                image.shape[1],
+                image.shape[0],
+                fill=False,
+                linestyle='dashed'))
+        ax.scatter(
+            mesh[:, 1],
+            mesh[:, 0],
+            c=values,
+            s=300,
+            cmap='gray',
+            edgecolors='r',
+            vmin=0,
+            vmax=1,
+            zorder=2)
+        for sp, mp in zip(sample_points, mesh):
+            ax.scatter(sp[:, 1], sp[:, 0], marker='x', color='blue', zorder=3)
+            if link_points:
+                for s in sp:
+                    ax.plot(
+                        [mp[1], s[1]], [mp[0], s[0]],
+                        color='red',
+                        alpha=.5,
+                        zorder=1)
+    return fig, axes
+
+
+if __name__ == '__main__':
+    ax = show_sampling()
+    plt.show()

nilearn/plotting/surf_plotting.py

@@ -6,15 +6,127 @@
 # Import libraries
 import nibabel
 import numpy as np
+
+# These will be removed (see comment on _points_in_unit_ball)
+from sklearn.externals import joblib
+from ..datasets.utils import _get_dataset_dir
+import sklearn.cluster
+# / These will be removed (see comment on _points_in_unit_ball)
+
 import matplotlib.pyplot as plt
 
 from mpl_toolkits.mplot3d import Axes3D
 from nibabel import gifti
 
+import nilearn.image
 from .._utils.compat import _basestring
+from .. import _utils
 from .img_plotting import _get_colorbar_and_data_ranges
 
 
+# Eventually, when this PR is ready, the sample locations inside the unit ball
+# will be hardcoded. For now we just compute them in a simple way (we will find
+# something better than the k-means hack) and cache the result.
+memory = joblib.Memory(_get_dataset_dir('joblib'), verbose=False)
+
+
+@memory.cache
+def _points_in_unit_ball(n_points=20, dim=3):
+    mc_cube = np.random.uniform(-1, 1, size=(5000, dim))
+    mc_ball = mc_cube[(mc_cube**2).sum(axis=1) <= 1.]
+    centroids, assignments, _ = sklearn.cluster.k_means(
+        mc_ball, n_clusters=n_points)
+    return centroids
+
+
+# Eventually will be replaced by nilearn.image.resampling.coord_transform, but
+# it only works for 3d images and 2d is useful for visu and debugging.
+def _transform_coord(coords, affine):
+    return affine.dot(np.vstack([coords.T, np.ones(coords.shape[0])]))[:-1].T
+
+
+def _ball_sample_locations(nodes, affine, ball_radius=3, n_points=20):
+    """Get n_points regularly spaced inside a ball."""
+    if affine is None:
+        affine = np.eye(nodes.shape[1] + 1)
+    offsets_world_space = _points_in_unit_ball(
+        dim=nodes.shape[1], n_points=n_points) * ball_radius
+    # later:
+    # offset_voxels = nilearn.image.resampling.coord_transform(
+    #     *offset_voxels.T, np.linalg.inv(affine))
+    mesh_voxel_space = _transform_coord(nodes, np.linalg.inv(affine))
+    linear_map = np.eye(affine.shape[0])
+    linear_map[:-1, :-1] = affine[:-1, :-1]
+    offsets_voxel_space = _transform_coord(offsets_world_space,
+                                           np.linalg.inv(linear_map))
+    sample_locations_voxel_space = (mesh_voxel_space[:, np.newaxis, :] +
+                                    offsets_voxel_space[np.newaxis, :])
+    return sample_locations_voxel_space, offsets_voxel_space
+
+
+def _ball_sampling(images, nodes, affine=None, ball_radius=3, n_points=20):
+    """In each image, average samples drawn from a ball around each node."""
+    images = np.asarray(images)
+    nodes = np.asarray(nodes)
+    sample_locations_voxel_space, offsets_voxel_space = _ball_sample_locations(
+        nodes, affine, ball_radius=ball_radius, n_points=n_points)
+    sample_indices = np.asarray(
+        np.floor(sample_locations_voxel_space), dtype=int)
+    pad_size = int(np.ceil(np.abs(offsets_voxel_space).max()))
+    pad = [(0, 0)] + [(pad_size, pad_size)] * nodes.shape[1]
+    sample_slices = [slice(None, None, None)] + np.s_[[
+        ax + pad_size for ax in np.rollaxis(sample_indices, -1)
+    ]]
+    try:
+        samples = np.pad(
+            images, pad, mode='constant',
+            constant_values=np.nan)[sample_slices]
+        texture = np.nanmean(samples, axis=2)
+        if not(np.isfinite(texture).all()):
+            raise IndexError()
+    except IndexError:
+        raise ValueError('Some nodes of the mesh are outside the image')
+    return texture, sample_indices, sample_locations_voxel_space
+
+
+def niimg_to_surf_data(image, mesh_nodes, ball_radius=3.):
+    """Extract surface data from a Nifti image.
+
+    Parameters
+    ----------
+
+    image : niimg-like object, 3d or 4d.
+
+    mesh_nodes : array-like,shape n_nodes * 3
+        The coordinates of the nodes of the mesh.
+
+    ball_radius : float, optional (default=3.).
+        The radius (in mm) of the ball over which image intensities are
+        averaged around each node.
+
+    Returns
+    -------
+    texture: array-like, 1d or 2d.
+        If image was a 3d image (e.g a stat map), a 1d vector is returned,
+        containing one value for each mesh node.
+        If image was a 4d image, a 2d array is returned, where each row
+        corresponds to a mesh node.
+
+    """
+    image = nilearn.image.load_img(image)
+    original_dimension = len(image.shape)
+    image = _utils.check_niimg(image, atleast_4d=True)
+    frames = np.rollaxis(image.get_data(), -1)
+    texture, indices, locations = _ball_sampling(
+        frames,
+        mesh_nodes,
+        _utils.compat.get_affine(image),
+        ball_radius=ball_radius)
+    if original_dimension == 3:
+        texture = texture[0]
+    return texture.T
+
+
 # function to figure out datatype and load data
 def load_surf_data(surf_data):
     """Loading data to be represented on a surface mesh.

nilearn/plotting/tests/test_surf_plotting.py

@@ -5,17 +5,24 @@
 
 from distutils.version import LooseVersion
 from nose import SkipTest
+from numpy.testing import (assert_array_equal, assert_array_almost_equal,
+                           assert_equal)
+from nose.tools import assert_true, assert_raises
+from nilearn._utils.testing import assert_raises_regex
 
-import matplotlib
 import numpy as np
-import nibabel as nb
+import matplotlib
 import matplotlib.pyplot as plt
+import sklearn.preprocessing
 
+import nibabel as nb
 from nibabel import gifti
-from numpy.testing import (assert_array_equal, assert_array_almost_equal,
-                           assert_equal)
 
-from nilearn._utils.testing import assert_raises_regex
+from nilearn import datasets
+import nilearn.image
+from nilearn.image.tests.test_resampling import rotation
+
+from nilearn.plotting import surf_plotting
 from nilearn.plotting.surf_plotting import (load_surf_data, load_surf_mesh,
                                             plot_surf, plot_surf_stat_map,
                                             plot_surf_roi)
@@ -388,3 +395,32 @@ def test_plot_surf_roi_error():
                         'Invalid input for roi_map',
                         plot_surf_roi, mesh,
                         roi_map={'roi1': roi1, 'roi2': roi2})
+
+
+def test_ball_sampling():
+    img = np.eye(3)
+    img[-1, -1] = 7.5
+    nodes = [[5, 5], [15, 5], [25, 25]]
+    affine = 10 * np.eye(3)
+    affine[-1, -1] = 1
+    texture = surf_plotting._ball_sampling(
+        [img], nodes, affine=affine, ball_radius=1)
+    assert_array_equal(texture[0][0], [1., 0., 7.5])
+    assert_raises(ValueError, surf_plotting._ball_sampling, [img], nodes)
+
+
+def test_niimg_to_surf_data():
+    mni = datasets.load_mni152_template()
+    fsaverage = datasets.fetch_surf_fsaverage5()
+    nodes = surf_plotting.load_surf_mesh(fsaverage.pial_left)[0]
+    proj_1 = surf_plotting.niimg_to_surf_data(mni, nodes)
+    assert_true(proj_1.ndim == 1)
+    mni_rot = nilearn.image.resample_img(
+        mni, target_affine=rotation(np.pi / 3., np.pi / 4.))
+    proj_2 = surf_plotting.niimg_to_surf_data(mni_rot, nodes)
+    assert_true((sklearn.preprocessing.normalize([proj_1])[0] *
+                 sklearn.preprocessing.normalize([proj_2])[0]).sum() > .998)
+    mni_4d = nilearn.image.concat_imgs([mni, mni])
+    proj_4d = surf_plotting.niimg_to_surf_data(mni_4d, nodes)
+    assert_array_equal(proj_4d.shape, [10242, 2])
+    assert_array_almost_equal(proj_4d[:, 0], proj_1, 3)