Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

LiFE model won't fit ? #2194

Closed
dPys opened this issue Jul 19, 2020 · 6 comments
Closed

LiFE model won't fit ? #2194

dPys opened this issue Jul 19, 2020 · 6 comments
Labels
type:Enhancement
Milestone
1.2

Comments

@dPys
Copy link
Contributor

dPys commented Jul 19, 2020 •
edited
Loading

@arokem -- Have you ever encountered the below error?

In [163]: fit_prob = fiber_model.fit(data, streamlines, affine=np.eye(4))                                                        
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-163-158d1438ad3f> in <module>
----> 1 fit_prob = fiber_model.fit(data, streamlines, affine=np.eye(4))

/usr/local/anaconda3/lib/python3.7/site-packages/dipy/tracking/life.py in fit(self, data, streamline, affine, evals, sphere)
    469             affine = np.eye(4)
    470         life_matrix, vox_coords = \
--> 471             self.setup(streamline, affine, evals=evals, sphere=sphere)
    472         (to_fit, weighted_signal, b0_signal, relative_signal, mean_sig,
    473          vox_data) = self._signals(data, vox_coords)

/usr/local/anaconda3/lib/python3.7/site-packages/dipy/tracking/life.py in setup(self, streamline, affine, evals, sphere)
    374         for s_idx, s in enumerate(streamline):
    375             if sphere is not False:
--> 376                 fiber_signal.append(SignalMaker.streamline_signal(s))
    377             else:
    378                 fiber_signal.append(streamline_signal(s, self.gtab, evals))

/usr/local/anaconda3/lib/python3.7/site-packages/dipy/tracking/life.py in streamline_signal(self, streamline)
    256         Approximate the signal for a given streamline
    257         """
--> 258         grad = streamline_gradients(streamline)
    259         sig_out = np.zeros((grad.shape[0], self.signal.shape[-1]))
    260         for ii, g in enumerate(grad):

/usr/local/anaconda3/lib/python3.7/site-packages/dipy/tracking/life.py in streamline_gradients(streamline)
    116 
    117     """
--> 118     return np.array(gradient(np.asarray(streamline))[0])
    119 
    120 

/usr/local/anaconda3/lib/python3.7/site-packages/dipy/tracking/life.py in gradient(f)
     81         slice3[axis] = 0
     82         # 1D equivalent -- out[0] = (f[1] - f[0])
---> 83         out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)])
     84         slice1[axis] = -1
     85         slice2[axis] = -1

IndexError: index 1 is out of bounds for axis 0 with size 1

The dwi and streamlines (depicted here as a density overlay) are aligned and the tractography is valid:
Screen Shot 2020-07-19 at 12 54 49 AM

Here's my gradient table:

In [171]: gtab.bvals                                                                                                             
Out[171]: 
array([   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.,    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.,    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.])

In [172]: gtab.bvecs                                                                                                             
Out[172]: 
array([[ 0.      ,  0.      ,  0.      ],
       [ 0.202244,  0.519251, -0.830347],
       [-0.199068,  0.519682, -0.830844],
       [-0.402972,  0.175034, -0.898319],
       [ 0.402691,  0.731728, -0.549922],
       [ 0.201685,  0.941542, -0.269855],
       [ 0.85303 ,  0.518601, -0.05825 ],
       [ 0.730986,  0.519821, -0.442092],
       [ 0.407118,  0.175802, -0.896297],
       [ 0.732969,  0.175604, -0.657206],
       [ 0.650673,  0.730677,  0.206729],
       [ 0.322091,  0.941382,  0.100288],
       [ 0.323097,  0.522558,  0.789013],
       [ 0.649862,  0.52144 ,  0.552974],
       [ 0.978606,  0.17679 ,  0.105244],
       [ 0.854464,  0.177024,  0.488419],
       [-0.001573,  0.733571,  0.679611],
       [-0.003094,  0.942801,  0.333342],
       [-0.655182,  0.52093 ,  0.547146],
       [-0.330365,  0.523133,  0.785615],
       [-0.196338, -0.179088, -0.964043],
       [-0.205773,  0.17956 ,  0.961985],
       [-0.651711,  0.731042,  0.202114],
       [-0.324485,  0.941187,  0.094218],
       [-0.199007,  0.940157, -0.27659 ],
       [-0.40215 ,  0.73038 , -0.552106],
       [-0.729851,  0.174471, -0.660967],
       [-0.728931,  0.516014, -0.449876],
       [-0.853626,  0.516529, -0.067227],
       [-0.859584,  0.175683,  0.479845],
       [-0.980035,  0.17419 ,  0.095863],
       [ 0.      ,  0.      ,  0.      ],
       [ 0.205649,  0.519283, -0.82949 ],
       [-0.195761,  0.52    , -0.831431],
       [-0.399543,  0.175277, -0.899802],
       [ 0.404589,  0.731911, -0.548284],
       [ 0.202122,  0.941763, -0.268753],
       [ 0.852773,  0.51948 , -0.054021],
       [ 0.732294,  0.521526, -0.437899],
       [ 0.411079,  0.177491, -0.894154],
       [ 0.735429,  0.17788 , -0.653837],
       [ 0.648873,  0.731173,  0.210596],
       [ 0.320369,  0.941574,  0.103932],
       [ 0.319563,  0.520786,  0.791619],
       [ 0.646901,  0.521262,  0.556601],
       [ 0.977867,  0.178346,  0.109404],
       [ 0.85218 ,  0.177812,  0.492109],
       [-0.005558,  0.731965,  0.681319],
       [-0.006584,  0.94195 ,  0.335689],
       [-0.658138,  0.51764 ,  0.546721],
       [-0.334001,  0.519749,  0.786323],
       [-0.19309 , -0.176899, -0.965102],
       [-0.208654,  0.176143,  0.961996],
       [-0.654037,  0.728732,  0.202943],
       [-0.327012,  0.940107,  0.096238],
       [-0.200403,  0.940519, -0.274339],
       [-0.402753,  0.731226, -0.550544],
       [-0.729872,  0.175417, -0.660693],
       [-0.728433,  0.51762 , -0.448838],
       [-0.853422,  0.516972, -0.066414],
       [-0.859794,  0.175172,  0.479654],
       [-0.980033,  0.174355,  0.09558 ],
       [ 0.      ,  0.      ,  0.      ],
       [ 0.206265,  0.52009 , -0.828832],
       [-0.195274,  0.520607, -0.831166],
       [-0.398807,  0.175695, -0.900047],
       [ 0.405043,  0.732445, -0.547233],
       [ 0.202429,  0.942155, -0.267147],
       [ 0.853092,  0.519035, -0.053252],
       [ 0.732637,  0.5213  , -0.437594],
       [ 0.411239,  0.17783 , -0.894013],
       [ 0.735455,  0.177448, -0.653925],
       [ 0.64948 ,  0.730624,  0.210629],
       [ 0.321488,  0.941179,  0.104054],
       [ 0.320455,  0.520275,  0.791595],
       [ 0.647746,  0.520589,  0.556249],
       [ 0.978147,  0.177104,  0.108912],
       [ 0.852591,  0.176851,  0.491744],
       [-0.004229,  0.732357,  0.680908],
       [-0.005135,  0.942253,  0.334861],
       [-0.657257,  0.519006,  0.546486],
       [-0.332945,  0.520814,  0.786066],
       [-0.19368 , -0.177309, -0.964909],
       [-0.208227,  0.177159,  0.961902],
       [-0.65334 ,  0.729499,  0.202427],
       [-0.325995,  0.94054 ,  0.095455],
       [-0.199654,  0.94055 , -0.274778],
       [-0.402167,  0.731248, -0.550943],
       [-0.729575,  0.175557, -0.660985],
       [-0.728305,  0.518085, -0.448508],
       [-0.853281,  0.517242, -0.066122],
       [-0.859912,  0.174809,  0.479576],
       [-0.980066,  0.174176,  0.095571]])

In [173]: gtab.b0s_mask                                                                                                          
Out[173]: 
array([ True, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False,  True, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False,  True,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False])
@arokem
Copy link
Contributor

arokem commented Jul 19, 2020 via email

@dPys
Copy link
Contributor Author

dPys commented Jul 20, 2020

@arokem -- Yep that appears to fixed it. Good call!
I've added a few additional precautions for input streamlines as well to ensure that 1) voxel indices are positive; 2) all streamline points are inside the B0 brain mask; and 3) minimum fiber length exceeds 10 mm.

I could also PR this if it'd be useful in Dipy, although it's mostly just an amalgamation of stuff that's already floating around? Let me know if you have any suggestions for further improvement:

def evaluate_streamline_plausibility(dwi_data, gtab, B0_mask_data, streamlines,
                                     affine=np.eye(4),
                                     sphere='repulsion724'):
    """
    Use Linear Fascicle Evaluation (LiFE) and other criteria to filter streamlines
    to only those with some evidence of neuroanatomical plausibility.

    Parameters
    ----------
    dwi_data : array
        4D array of dwi data.
    gtab : Obj
        DiPy object storing diffusion gradient information.
    B0_mask_data : array
       B0 brain mask 3d data.
    streamlines : ArraySequence
        DiPy list/array-like object of streamline points from tractography.

    Returns
    -------
    streamlines : ArraySequence
        DiPy list/array-like object of streamline points from tractography.

    References
    ----------
    .. [1] Pestilli, F., Yeatman, J, Rokem, A. Kay, K. and Wandell B.A. (2014).
     Validation and statistical inference in living connectomes.
     Nature Methods 11: 1058-1063. doi:10.1038/nmeth.3098
    """
    import dipy.tracking.life as life
    import dipy.core.optimize as opt
    from dipy.tracking._utils import _mapping_to_voxel
    from dipy.data import get_sphere
    from dipy.tracking import utils
    from dipy.tracking.streamline import Streamlines

    print('Filtering streamlines by length > 10...')
    streamlines_long = nib.streamlines. \
        array_sequence.ArraySequence(
        [
            s
            for s in streamlines
            if len(s) >= float(10)
        ]
    )
    print('Removing streamlines with negative voxel indices...')
    # Remove any streamlines with negative voxel indices
    lin_T, offset = _mapping_to_voxel(np.eye(4))
    streamlines_positive = []
    for sl in streamlines_long:
        inds = np.dot(sl, lin_T)
        inds += offset
        if not inds.min().round(decimals=6) < 0:
            streamlines_positive.append(sl)
    del streamlines_long
    # Filter resulting streamlines by those that stay entirely
    # inside the brain
    streamlines_in_brain = Streamlines(utils.target(
        streamlines_positive, np.eye(4),
        B0_mask_data.astype('bool'), include=True
    ))
    streamlines_in_brain = [i for i in streamlines_in_brain]
    del streamlines_positive
    print('Performing Linear Fascicle Evaluation...')
    original_count = len(streamlines)
    sphere = get_sphere(sphere)
    fiber_model = life.FiberModel(gtab)
    fiber_fit = fiber_model.fit(dwi_data, streamlines_in_brain,
                                affine=affine,
                                sphere=sphere)
    streamlines = list(np.array(streamlines_in_brain)[
                           np.where(fiber_fit.beta > 0)[0]])
    pruned_count = len(streamlines)
    model_predict = fiber_fit.predict()
    model_error = model_predict - fiber_fit.data
    model_rmse = np.sqrt(np.mean(model_error[:, 10:] ** 2, -1))
    beta_baseline = np.zeros(fiber_fit.beta.shape[0])
    pred_weighted = np.reshape(opt.spdot(fiber_fit.life_matrix, beta_baseline),
                               (fiber_fit.vox_coords.shape[0],
                                np.sum(~gtab.b0s_mask)))
    mean_pred = np.empty((fiber_fit.vox_coords.shape[0], gtab.bvals.shape[0]))
    S0 = fiber_fit.b0_signal
    mean_pred[..., gtab.b0s_mask] = S0[:, None]
    mean_pred[..., ~gtab.b0s_mask] = (pred_weighted +
                                      fiber_fit.mean_signal
                                      [:, None]) * S0[:, None]
    mean_error = mean_pred - fiber_fit.data
    mean_rmse = np.sqrt(np.mean(mean_error ** 2, -1))
    print(f"Original # Streamlines: {original_count}")
    print(f"Final # Streamlines: {pruned_count}")
    print(f"Streamlines removed: {pruned_count - original_count}")
    print(f"Mean RMSE: {np.mean(mean_rmse)}")
    print(f"Mean Model RMSE: {np.mean(model_rmse)}")
    print(f"Mean Reduction RMSE: {np.mean(mean_rmse - model_rmse)}")
    return streamlines

@dPys
Copy link
Contributor Author

dPys commented Jul 20, 2020

@Garyfallidis It would still be nice if there was a way to incorporate further plausibility check on the basis of known bundles? but I haven't been able to find a good solution that is fast enough and doesn't require some degree of tuning depending on which fibers in particular are being tracked. Could it be feasible to maybe even just roughly evaluate similarity between a density map of the streamlines and, say, a probabilistic white-matter atlas resampled into the same space?

@arokem
Copy link
Contributor

arokem commented Jul 20, 2020

For starters, could you add just some simple input validation in the fit method, making sure that all the streamlines are longer than 1 node?

@dPys
Copy link
Contributor Author

dPys commented Jul 20, 2020

#2195 @arokem

@skoudoro
Copy link
Member

skoudoro commented Sep 1, 2020

I supposed it is fixed by #2195. Closing

@skoudoro skoudoro closed this as completed Sep 1, 2020
@skoudoro skoudoro added this to the 1.2 milestone Sep 1, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
type:Enhancement
Projects
None yet
Milestone
1.2
Development

No branches or pull requests
3 participants
@arokem @dPys @skoudoro