1
2import os
3import numpy as np
4from scipy.spatial import cKDTree
5from scipy.ndimage.interpolation import map_coordinates
6from scipy.spatial.distance import mahalanobis
7
8from dipy.utils.optpkg import optional_package
9from dipy.io.utils import save_buan_profiles_hdf5
10from dipy.segment.clustering import QuickBundles
11from dipy.segment.metric import AveragePointwiseEuclideanMetric
12from dipy.tracking.streamline import (set_number_of_points,
13                                      values_from_volume,
14                                      orient_by_streamline,
15                                      transform_streamlines,
16                                      Streamlines)
17
18pd, have_pd, _ = optional_package("pandas")
19_, have_tables, _ = optional_package("tables")
20
21
22def peak_values(bundle, peaks, dt, pname, bname, subject, group_id, ind, dir):
23    """ Peak_values function finds the generalized fractional anisotropy (gfa)
24        and quantitative anisotropy (qa) values from peaks object (eg: csa) for
25        every point on a streamline used while tracking and saves it in hd5
26        file.
27
28    Parameters
29    ----------
30    bundle : string
31        Name of bundle being analyzed
32    peaks : peaks
33        contains peak directions and values
34    dt : DataFrame
35        DataFrame to be populated
36    pname : string
37        Name of the dti metric
38    bname : string
39        Name of bundle being analyzed.
40    subject : string
41        subject number as a string (e.g. 10001)
42    group_id : integer
43        which group subject belongs to 1 patient and 0 for control
44    ind : integer list
45        ind tells which disk number a point belong.
46    dir : string
47        path of output directory
48
49    """
50
51    gfa = peaks.gfa
52    anatomical_measures(bundle, gfa, dt, pname+'_gfa', bname, subject, group_id,
53                        ind, dir)
54
55    qa = peaks.qa[...,0]
56    anatomical_measures(bundle, qa, dt, pname+'_qa', bname, subject, group_id,
57                        ind, dir)
58
59
60def anatomical_measures(bundle, metric, dt, pname, bname, subject, group_id,
61                        ind, dir):
62    """ Calculates dti measure (eg: FA, MD) per point on streamlines and
63        save it in hd5 file.
64
65    Parameters
66    ----------
67    bundle : string
68        Name of bundle being analyzed
69    metric : matrix of float values
70        dti metric e.g. FA, MD
71    dt : DataFrame
72        DataFrame to be populated
73    pname : string
74        Name of the dti metric
75    bname : string
76        Name of bundle being analyzed.
77    subject : string
78        subject number as a string (e.g. 10001)
79    group_id : integer
80        which group subject belongs to 1 for patient and 0 control
81    ind : integer list
82        ind tells which disk number a point belong.
83    dir : string
84        path of output directory
85
86    """
87    dt["streamline"] = []
88    dt["disk"] = []
89    dt["subject"] = []
90    dt[pname] = []
91    dt["group"] = []
92
93    values = map_coordinates(metric, bundle._data.T,
94                             order=1)
95
96    dt["disk"].extend(ind[list(range(len(values)))]+1)
97    dt["subject"].extend([subject]*len(values))
98    dt["group"].extend([group_id]*len(values))
99    dt[pname].extend(values)
100
101    for st_i in range(len(bundle)):
102
103        st = bundle[st_i]
104        dt["streamline"].extend([st_i]*len(st))
105
106    file_name = bname + "_" + pname
107
108    save_buan_profiles_hdf5(os.path.join(dir, file_name), dt)
109
110
111def assignment_map(target_bundle, model_bundle, no_disks):
112    """
113    Calculates assignment maps of the target bundle with reference to
114    model bundle centroids.
115
116    Parameters
117    ----------
118    target_bundle : streamlines
119        target bundle extracted from subject data in common space
120    model_bundle : streamlines
121        atlas bundle used as reference
122    no_disks : integer, optional
123        Number of disks used for dividing bundle into disks. (Default 100)
124
125    References
126    ----------
127    .. [Chandio2020] Chandio, B.Q., Risacher, S.L., Pestilli, F., Bullock, D.,
128    Yeh, FC., Koudoro, S., Rokem, A., Harezlak, J., and Garyfallidis, E.
129    Bundle analytics, a computational framework for investigating the
130    shapes and profiles of brain pathways across populations.
131    Sci Rep 10, 17149 (2020)
132
133    """
134
135    mbundle_streamlines = set_number_of_points(model_bundle,
136                                               nb_points=no_disks)
137
138    metric = AveragePointwiseEuclideanMetric()
139    qb = QuickBundles(threshold=85., metric=metric)
140    clusters = qb.cluster(mbundle_streamlines)
141    centroids = Streamlines(clusters.centroids)
142
143    _, indx = cKDTree(centroids.get_data(), 1,
144                      copy_data=True).query(target_bundle.get_data(), k=1)
145
146    return indx
147
148
149def gaussian_weights(bundle, n_points=100, return_mahalnobis=False,
150                     stat=np.mean):
151    """
152    Calculate weights for each streamline/node in a bundle, based on a
153    Mahalanobis distance from the core the bundle, at that node (mean, per
154    default).
155
156    Parameters
157    ----------
158    bundle : Streamlines
159        The streamlines to weight.
160    n_points : int, optional
161        The number of points to resample to. *If the `bundle` is an array, this
162        input is ignored*. Default: 100.
163
164    Returns
165    -------
166    w : array of shape (n_streamlines, n_points)
167        Weights for each node in each streamline, calculated as its relative
168        inverse of the Mahalanobis distance, relative to the distribution of
169        coordinates at that node position across streamlines.
170
171    """
172    # Resample to same length for each streamline:
173    bundle = set_number_of_points(bundle, n_points)
174
175    # This is the output
176    w = np.zeros((len(bundle), n_points))
177
178    # If there's only one fiber here, it gets the entire weighting:
179    if len(bundle) == 1:
180        if return_mahalnobis:
181            return np.array([np.nan])
182        else:
183            return np.array([1])
184
185    for node in range(n_points):
186        # This should come back as a 3D covariance matrix with the spatial
187        # variance covariance of this node across the different streamlines
188        # This is a 3-by-3 array:
189        node_coords = bundle._data[node::n_points]
190        c = np.cov(node_coords.T, ddof=0)
191        # Reorganize as an upper diagonal matrix for expected Mahalanobis
192        # input:
193        c = np.array([[c[0, 0], c[0, 1], c[0, 2]],
194                      [0, c[1, 1], c[1, 2]],
195                      [0, 0, c[2, 2]]])
196        # Calculate the mean or median of this node as well
197        # delta = node_coords - np.mean(node_coords, 0)
198        m = stat(node_coords, 0)
199        # Weights are the inverse of the Mahalanobis distance
200        for fn in range(len(bundle)):
201            # In the special case where all the streamlines have the exact same
202            # coordinate in this node, the covariance matrix is all zeros, so
203            # we can't calculate the Mahalanobis distance, we will instead give
204            # each streamline an identical weight, equal to the number of
205            # streamlines:
206            if np.allclose(c, 0):
207                w[:, node] = len(bundle)
208                break
209            # Otherwise, go ahead and calculate Mahalanobis for node on
210            # fiber[fn]:
211            w[fn, node] = mahalanobis(node_coords[fn], m, np.linalg.inv(c))
212    if return_mahalnobis:
213        return w
214    # weighting is inverse to the distance (the further you are, the less you
215    # should be weighted)
216    w = 1 / w
217    # Normalize before returning, so that the weights in each node sum to 1:
218    return w / np.sum(w, 0)
219
220
221def afq_profile(data, bundle, affine, n_points=100, profile_stat=np.average,
222                orient_by=None, weights=None, **weights_kwarg):
223    """
224    Calculates a summarized profile of data for a bundle or tract
225    along its length.
226
227    Follows the approach outlined in [Yeatman2012]_.
228
229    Parameters
230    ----------
231    data : 3D volume
232        The statistic to sample with the streamlines.
233
234    bundle : StreamLines class instance
235        The collection of streamlines (possibly already resampled into an array
236         for each to have the same length) with which we are resampling. See
237         Note below about orienting the streamlines.
238    affine : array_like (4, 4)
239        The mapping from voxel coordinates to streamline points.
240        The voxel_to_rasmm matrix, typically from a NIFTI file.
241    n_points: int, optional
242        The number of points to sample along the bundle. Default: 100.
243    orient_by: streamline, optional.
244        A streamline to use as a standard to orient all of the streamlines in
245        the bundle according to.
246    weights : 1D array or 2D array or callable (optional)
247        Weight each streamline (1D) or each node (2D) when calculating the
248        tract-profiles. Must sum to 1 across streamlines (in each node if
249        relevant). If callable, this is a function that calculates weights.
250    profile_stat : callable
251        The statistic used to average the profile across streamlines.
252        If weights is not None, this must take weights as a keyword argument.
253        The default, np.average, is the same as np.mean but takes weights
254        as a keyword argument.
255    weights_kwarg : key-word arguments
256        Additional key-word arguments to pass to the weight-calculating
257        function. Only to be used if weights is a callable.
258
259    Returns
260    -------
261    ndarray : a 1D array with the profile of `data` along the length of
262        `bundle`
263
264    Notes
265    -----
266    Before providing a bundle as input to this function, you will need to make
267    sure that the streamlines in the bundle are all oriented in the same
268    orientation relative to the bundle (use :func:`orient_by_streamline`).
269
270    References
271    ----------
272    .. [Yeatman2012] Yeatman, Jason D., Robert F. Dougherty,
273       Nathaniel J. Myall, Brian A. Wandell, and Heidi M. Feldman. 2012.
274       "Tract Profiles of White Matter Properties: Automating Fiber-Tract
275       Quantification" PloS One 7 (11): e49790.
276
277    """
278    if orient_by is not None:
279        bundle = orient_by_streamline(bundle, orient_by)
280    if affine is None:
281        affine = np.eye(4)
282    if len(bundle) == 0:
283        raise ValueError("The bundle contains no streamlines")
284
285    # Resample each streamline to the same number of points:
286    fgarray = set_number_of_points(bundle, n_points)
287
288    # Extract the values
289    values = np.array(values_from_volume(data, fgarray, affine))
290
291    if weights is not None:
292        if callable(weights):
293            weights = weights(bundle, **weights_kwarg)
294        else:
295            # We check that weights *always sum to 1 across streamlines*:
296            if not np.allclose(np.sum(weights, 0), np.ones(n_points)):
297                raise ValueError(
298                    "The sum of weights across streamlines",
299                    " must be equal to 1")
300
301        return profile_stat(values, weights=weights, axis=0)
302    else:
303        return profile_stat(values, axis=0)
304