Merge pull request #96 from jhlegarreta/AddGradientPlotMethod

ENH: Add gradient plot method

nireports/conftest.py

@@ -33,7 +33,7 @@
 # disable ET
 os.environ['NO_ET'] = '1'
 
-_datadir = (Path(__file__).parent / "tests" / "data").resolve(strict=True)
+_datadir = (Path(__file__).parent / "tests" / "data").absolute()
 niprepsdev_path = os.getenv(
     "TEST_DATA_HOME", str(Path.home() / ".cache" / "nipreps-dev")
 )

nireports/reportlets/modality/dwi.py

@@ -23,6 +23,8 @@
 """Visualizations for diffusion MRI data."""
 import numpy as np
 from matplotlib import pyplot as plt
+from matplotlib.pyplot import cm
+from mpl_toolkits.mplot3d import art3d
 
 
 def plot_heatmap(
@@ -124,3 +126,221 @@ def plot_heatmap(
     fig.tight_layout(rect=[0.02, 0, 1, 1])
 
     return fig
+
+
+
+def rotation_matrix(u, v):
+    r"""
+    Calculate the rotation matrix *R* such that :math:`R \cdot \mathbf{u} = \mathbf{v}`.
+
+    Extracted from `Emmanuel Caruyer's code
+    <https://github.com/ecaruyer/qspace/blob/master/qspace/visu/visu_points.py>`__,
+    which is distributed under the revised BSD License:
+    Copyright (c) 2013-2015, Emmanuel Caruyer
+    All rights reserved.
+
+    .. admonition :: List of changes
+        Only minimal updates to leverage Numpy.
+
+    Parameters
+    ----------
+    u : :obj:`numpy.ndarray`
+        A vector.
+    v : :obj:`numpy.ndarray`
+        A vector.
+
+    Returns
+    -------
+    R : :obj:`numpy.ndarray`
+        The rotation matrix.
+
+    """
+
+    # the axis is given by the product u x v
+    u = u / np.linalg.norm(u)
+    v = v / np.linalg.norm(v)
+    w = np.asarray(
+        [
+            u[1] * v[2] - u[2] * v[1],
+            u[2] * v[0] - u[0] * v[2],
+            u[0] * v[1] - u[1] * v[0],
+        ]
+    )
+    if (w ** 2).sum() < (np.finfo(w.dtype).eps * 10):
+        # The vectors u and v are collinear
+        return np.eye(3)
+
+    # Compute sine and cosine
+    c = u @ v
+    s = np.linalg.norm(w)
+
+    w = w / s
+    P = np.outer(w, w)
+    Q = np.asarray([[0, -w[2], w[1]], [w[2], 0, -w[0]], [-w[1], w[0], 0]])
+    R = P + c * (np.eye(3) - P) + s * Q
+    return R
+
+
+def draw_circles(positions, radius, n_samples=20):
+    r"""
+    Draw circular patches (lying on a sphere) at given positions.
+
+    Adapted from `Emmanuel Caruyer's code
+    <https://github.com/ecaruyer/qspace/blob/master/qspace/visu/visu_points.py>`__,
+    which is distributed under the revised BSD License:
+    Copyright (c) 2013-2015, Emmanuel Caruyer
+    All rights reserved.
+
+    .. admonition :: List of changes
+        Modified to take the full list of normalized bvecs and corresponding circle
+        radii instead of taking the list of bvecs and radii for a specific shell
+        (*b*-value).
+
+    Parameters
+    ----------
+    positions : :obj:`numpy.ndarray`
+        An array :math:`N \times 3` of 3D cartesian positions.
+    radius : :obj:`float`
+        The reference radius (or, the radius in single-shell plots)
+    n_samples : :obj:`int`
+        The number of samples on the sphere.
+
+    Returns
+    -------
+    circles : :obj:`numpy.ndarray`
+        Circular patches.
+
+    """
+
+    # A circle centered at [1, 0, 0] with radius r
+    t = np.linspace(0, 2 * np.pi, n_samples)
+
+    nb_points = positions.shape[0]
+    circles = np.zeros((nb_points, n_samples, 3))
+    for i in range(positions.shape[0]):
+        circle_x = np.zeros((n_samples, 3))
+        dots_radius = np.sqrt(radius[i]) * 0.04
+        circle_x[:, 1] = dots_radius * np.cos(t)
+        circle_x[:, 2] = dots_radius * np.sin(t)
+        norm = np.linalg.norm(positions[i])
+        point = positions[i] / norm
+        r1 = rotation_matrix(np.asarray([1, 0, 0]), point)
+        circles[i] = positions[i] + np.dot(r1, circle_x.T).T
+    return circles
+
+
+def draw_points(gradients, ax, rad_min=0.3, rad_max=0.7, colormap="viridis"):
+    """
+    Draw the vectors on a shell.
+
+    Adapted from `Emmanuel Caruyer's code
+    <https://github.com/ecaruyer/qspace/blob/master/qspace/visu/visu_points.py>`__,
+    which is distributed under the revised BSD License:
+    Copyright (c) 2013-2015, Emmanuel Caruyer
+    All rights reserved.
+
+    .. admonition :: List of changes
+        * The input is a single 2D numpy array of the gradient table in RAS+B format
+        * The scaling of the circle radius for each bvec proportional to the inverse of
+          the bvals. A minimum/maximal value for the radii can be specified.
+        * Circles for each bvec are drawn at once instead of looping over the shells.
+        * Some variables have been renamed (like vects to bvecs)
+
+    Parameters
+    ----------
+    gradients : :obj:`numpy.ndarray`
+        An (N, 4) shaped array of the gradient table in RAS+B format.
+    ax : :obj:`matplotlib.axes.Axis`
+        The matplotlib axes instance to plot in.
+    rad_min : :obj:`float` between 0 and 1
+        Minimum radius of the circle that renders a gradient direction.
+    rad_max : :obj:`float` between 0 and 1
+        Maximum radius of the circle that renders a gradient direction.
+    colormap : :obj:`matplotlib.pyplot.cm.ColorMap`
+        matplotlib colormap name.
+
+    """
+
+    # Initialize 3D view
+    elev = 90
+    azim = 0
+    ax.view_init(azim=azim, elev=elev)
+
+    # Normalize to 1 the highest bvalue
+    bvals = np.copy(gradients[:, 3])
+    bvals = bvals / bvals.max()
+
+    # Colormap depending on bvalue (for visualization)
+    cmap = cm.get_cmap(colormap)
+    colors = cmap(bvals)
+
+    # Relative shell radii proportional to the inverse of bvalue (for visualization)
+    rs = np.reciprocal(bvals)
+    rs = rs / rs.max()
+
+    # Readjust radius of the circle given the minimum and maximal allowed values.
+    if rs.min() != rs.max():
+        rs = rs - rs.min()
+        rs = rs / (rs.max() - rs.min())
+    rs = rs * (rad_max - rad_min) + rad_min
+
+    bvecs = np.copy(
+        gradients[:, :3],
+    )
+    bvecs[bvecs[:, 2] < 0] *= -1
+
+    # Render all gradient direction of all b-values
+    circles = draw_circles(bvecs, rs)
+    ax.add_collection(art3d.Poly3DCollection(circles, facecolors=colors, linewidth=0))
+
+    max_val = 0.6
+    ax.set_xlim(-max_val, max_val)
+    ax.set_ylim(-max_val, max_val)
+    ax.set_zlim(-max_val, max_val)
+    ax.axis("off")
+
+
+def plot_gradients(
+    gradients,
+    title=None,
+    ax=None,
+    spacing=0.05,
+    **kwargs,
+):
+    """
+    Draw the vectors on a unit sphere with color code for multiple b-value.
+
+    Parameters
+    ----------
+    gradients : :obj:`numpy.ndarray`
+        An (N, 4) shaped array of the gradient table in RAS+B format.
+    title : :obj:`str`
+        Plot title.
+    ax : :obj:`matplotlib.axes.Axis`
+        A figure's axis to plot on.
+    spacing : :obj:`float`
+        Plot spacing.
+    kwargs : :obj:`dict`
+        Extra args given to :obj:`eddymotion.viz.draw_points()`.
+
+    Returns
+    -------
+    ax : :obj:`matplotlib.axes.Axis`
+        The figure's axis where the data is plotted.
+
+    """
+
+    # Initialize figure
+    if ax is None:
+        figsize = kwargs.pop("figsize", (9.0, 9.0))
+        fig = plt.figure(figsize=figsize)
+        ax = fig.add_subplot(111, projection="3d")
+    plt.subplots_adjust(bottom=spacing, top=1 - spacing, wspace=2 * spacing)
+
+    # Draw points after re-projecting all shells to the unit sphere
+    draw_points(gradients, ax, **kwargs)
+
+    if title:
+        plt.suptitle(title)
+
+    return ax

nireports/tests/data/ds000114_singleshell.bval

@@ -0,0 +1 @@
+0 0 0 0 0 0 0 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 

nireports/tests/data/ds000114_singleshell.bvec

@@ -0,0 +1,3 @@
+0 0 0 0 0 0 0 -1 -0.002 0.026 -0.591 0.236 0.893 -0.796 -0.234 -0.936 -0.506 -0.346 -0.457 0.487 0.618 0.577 0.827 -0.894 -0.29 -0.116 0.8 -0.514 0.789 -0.949 -0.233 0.021 -0.217 -0.774 0.161 0.147 -0.888 0.562 0.381 0.306 0.332 0.963 0.959 -0.453 0.773 -0.709 0.693 -0.682 0.142 0.74 0.103 -0.584 0.088 0.552 -0.838 -0.363 0.184 0.721 -0.433 -0.502 0.171 -0.463 -0.385 0.713 -0.26 -0.001 -0.037 -0.57 0.282 -0.721 -0.267 
+0 0 0 0 0 0 0 0 1 0.649 -0.766 -0.524 -0.259 0.129 0.93 0.14 -0.845 -0.847 -0.631 -0.389 0.673 -0.105 -0.521 -0.04 -0.541 -0.963 0.403 0.84 0.153 -0.233 0.783 -0.188 -0.956 -0.604 0.356 0.731 0.417 0.232 0.143 -0.199 -0.13 -0.265 0.205 -0.889 0.628 0.408 0.024 0.529 -0.725 0.388 0.822 -0.596 -0.335 -0.792 -0.458 -0.561 0.392 -0.693 0.682 0.69 -0.509 0.423 -0.809 -0.247 0.885 0.077 -0.902 -0.303 0.145 0.608 0.96 
+0 0 0 0 0 0 0 0 0 0.76 0.252 0.818 0.368 0.591 0.284 0.324 -0.175 -0.402 -0.627 0.782 0.407 -0.81 0.213 -0.447 -0.789 -0.245 -0.444 0.174 -0.596 0.211 0.577 -0.982 0.199 0.19 0.921 -0.666 0.193 -0.794 0.914 -0.931 0.934 0.044 0.193 0.068 0.088 0.575 0.721 -0.506 0.674 0.549 0.56 0.551 0.938 0.259 -0.296 0.744 -0.901 0.009 -0.589 0.521 -0.844 0.779 0.444 0.656 -0.387 -0.997 0.43 -0.763 -0.948 0.332 -0.085 

nireports/tests/data/ds004737_dsi.bval

@@ -0,0 +1 @@
+5 5 4195 3590 4790 4990 3790 2590 1800 4195 995 4990 2595 2795 5 5000 3395 1795 3595 2790 2190 4390 5000 1195 1195 2795 4990 5 1795 995 2195 4795 4985 5000 2795 1790 1595 4795 1995 800 5 2790 2190 3790 5000 1595 4990 2795 1195 2195 1800 2000 3595 5 4395 2000 2600 4990 2000 2195 4390 800 2795 1595 2790 4790 5 3395 1195 1795 2195 5000 4990 2395 4795 4795 1000 1195 995 5 2190 1990 4990 2795 3395 995 5000 2595 3400 1990 795 995 1795 1990 1195 1795 1795 2600 2395 4990 1000 1600 1995 5

nireports/tests/data/ds004737_dsi.bvec

@@ -0,0 +1,3 @@
+-0.5 -0.5 -0.873381 0.707894 0.408703 0 0.229712 0 0 -0.436552 0.895769 -0.800629 0 0.802493 -0.5 -0.599999 -0.970507 0.333952 -0.943136 -0.267753 -0.302149 0.426759 0.999999 -0.408871 0.817403 0.534767 0 -0.5 0.333952 0 0.301735 0.408432 0 -0.799999 -0.802199 0 -0.708215 0.408432 0 0 -0.5 -0.53522 -0.302149 -0.229712 0 0 -0.600667 -0.534767 0.817403 -0.301735 0.999998 -0.948681 -0.943136 -0.5 -0.639981 0.948681 0.832049 0.600667 0.316227 -0.301735 -0.426759 0.999993 0.802199 0.708215 0.267753 0.408703 -0.5 0.485415 -0.409158 -0.667533 -0.905069 0.799999 0 -0.577993 -0.408432 -0.81695 0.894423 0.408871 0.448471 -0.5 0.302149 0 0.800629 0.534767 0.970507 -0.895769 0.599999 0 0.242535 -0.31699 0 0 0.667533 0 -0.409158 -0.667533 -0.333952 -0.554699 -0.577993 0 -0.894423 -0.707104 0.94933 -0.5
+-0.5 -0.5 0.218371 0 0.408703 -0.600667 -0.688954 0.555515 -0.999998 -0.873155 1e-08 0 -0.832872 0.267561 -0.5 -0.799999 0 0.667533 -0.2358 0.53522 0.302149 -0.640211 0 0.817403 -0.408871 -0.802199 -0.800629 -0.5 -0.667533 0.448471 0.905069 0.81695 0 -0.6 0.534767 0 0 -0.81695 0.94933 -0.999993 -0.5 -0.267753 -0.302149 -0.688954 -0.999999 -0.708215 0 -0.802199 0.408871 -0.905069 0 -0.316227 0.2358 -0.5 0.639981 -0.316227 -0.554699 0 -0.948681 0.905069 0.640211 0 0.534767 0 0.53522 -0.408703 -0.5 -0.728195 0.409158 0.333952 0.301735 -0.6 0.800629 0.577993 0.81695 -0.408433 -0.44721 -0.817403 0 -0.5 -0.302149 0.31699 0 0.802199 0 -1e-08 -0.799999 0.832872 -0.970142 0 0 -0.895769 0.333952 -0.31699 -0.409158 -0.333952 -0.667533 -0.832049 -0.577993 0.600667 -0.44721 -0.707104 0 -0.5
+0.707107 0.707107 -0.435338 -0.706319 -0.816042 -0.799499 -0.687441 -0.831507 0.00222221 -0.216846 -0.44452 -0.59916 -0.553466 -0.533306 0.707107 0.00128 -0.241074 -0.665489 -0.234294 -0.801154 -0.904109 -0.638754 0.0014 -0.405803 -0.405803 -0.265521 -0.59916 0.707107 -0.665489 -0.893797 -0.299677 -0.407155 -1 0.00128 -0.265521 -1 -0.705997 -0.407155 -0.31428 0.00374992 0.707107 -0.801154 -0.904109 -0.687441 0.0014 -0.705997 -0.799499 -0.265521 -0.405803 -0.299677 0.00222221 0.00205547 -0.234294 0.707107 -0.425263 0.00205547 0.00192011 -0.799499 0.00205547 -0.299677 -0.638754 0.00374992 -0.265521 -0.705997 -0.801154 -0.816042 0.707107 -0.483843 -0.815585 -0.665489 -0.299677 0.00128 -0.59916 -0.576063 -0.407155 -0.407155 0.00313045 -0.405803 -0.893797 0.707107 -0.904109 -0.948429 -0.59916 -0.265521 -0.241074 -0.44452 0.00128 -0.553466 0.00149801 -0.948429 -1 -0.44452 -0.665489 -0.948429 -0.815585 -0.665489 -0.665489 0.00192011 -0.576063 -0.799499 0.00313045 0.00265162 -0.31428 0.707107

nireports/tests/data/hcph_multishell.bval

@@ -0,0 +1 @@
+0 3000 3000 3000 2000 2000 700 3000 2000 1000 1000 2000 3000 3000 1000 1000 3000 3000 3000 3000 2000 3000 2000 2000 3000 2000 2000 2000 1000 3000 3000 3000 3000 3000 3000 2000 3000 2000 2000 2000 3000 2000 3000 2000 1000 3000 2000 0 3000 1000 3000 3000 2000 2000 700 3000 3000 2000 3000 700 3000 2000 2000 3000 2000 2000 3000 2000 3000 3000 2000 3000 2000 3000 1000 1000 3000 1000 2000 700 2000 1000 2000 3000 3000 700 3000 2000 3000 1000 3000 2000 1000 2000 0 3000 1000 1000 3000 3000 3000 1000 3000 2000 2000 2000 3000 3000 2000 2000 3000 700 3000 3000 2000 1000 2000 3000 3000 3000 2000 3000 3000 3000 1000 3000 2000 1000 2000 2000 2000 3000 3000 3000 1000 3000 3000 3000 3000 3000 0 3000 3000 3000 2000 2000 700 3000 2000 1000 1000 2000 3000 3000 1000 1000 3000 3000 3000 3000 2000 3000 2000 2000 3000 2000 2000 2000 1000 3000 3000 3000 3000 3000 3000 2000 3000 2000 2000 2000 3000 2000 3000 2000 1000 3000 2000 0 3000 1000 3000 3000 2000 2000 700 3000 3000 2000 3000 700 3000 2000 2000 3000 2000 2000 3000 2000 3000 3000 2000 3000 2000 3000 1000 1000 3000 1000 2000 700 2000 1000 2000 3000 3000 700 3000 2000 3000 1000 3000 2000 1000 2000 0 3000 1000 1000 3000 3000 3000 1000 3000 2000 2000 2000 3000 3000 2000 2000 3000 700 3000 3000 2000 1000 2000 3000 3000 3000 2000 3000 3000 3000 1000 3000 2000 1000 2000 2000 2000 3000 3000 3000 1000 3000 3000 3000 3000 3000