msd response function

fs_keys.p

(dp1
S'cerebellum'
p2
cnumpy.core.multiarray
_reconstruct
p3
(cnumpy
ndarray
p4
(I0
tS'b'
tRp5
(I1
(I4
tcnumpy
dtype
p6
(S'i8'
I0
I1
tRp7
(I3
S'<'
NNNI-1
I-1
I0
tbI00
S'\x07\x00\x00\x00\x00\x00\x00\x00\x08\x00\x00\x00\x00\x00\x00\x00.\x00\x00\x00\x00\x00\x00\x00/\x00\x00\x00\x00\x00\x00\x00'
tbsS'Unknown'
p8
g3
(g4
(I0
tS'b'
tRp9
(I1
(I5
tg7
I00
S'\x00\x00\x00\x00\x00\x00\x00\x00M\x00\x00\x00\x00\x00\x00\x00P\x00\x00\x00\x00\x00\x00\x00\xe8\x03\x00\x00\x00\x00\x00\x00\xd0\x07\x00\x00\x00\x00\x00\x00'
tbsS'wm'
p10
g3
(g4
(I0
tS'b'
tRp11
(I1
(I8
tg7
I00
S'\x02\x00\x00\x00\x00\x00\x00\x00)\x00\x00\x00\x00\x00\x00\x00U\x00\x00\x00\x00\x00\x00\x00\xfb\x00\x00\x00\x00\x00\x00\x00\xfc\x00\x00\x00\x00\x00\x00\x00\xfd\x00\x00\x00\x00\x00\x00\x00\xfe\x00\x00\x00\x00\x00\x00\x00\xff\x00\x00\x00\x00\x00\x00\x00'
tbsS'csf'
p12
g3
(g4
(I0
tS'b'
tRp13
(I1
(I12
tg7
I00
S'\x04\x00\x00\x00\x00\x00\x00\x00\x05\x00\x00\x00\x00\x00\x00\x00\x0e\x00\x00\x00\x00\x00\x00\x00\x0f\x00\x00\x00\x00\x00\x00\x00\x18\x00\x00\x00\x00\x00\x00\x00\x1e\x00\x00\x00\x00\x00\x00\x00\x1f\x00\x00\x00\x00\x00\x00\x00+\x00\x00\x00\x00\x00\x00\x00,\x00\x00\x00\x00\x00\x00\x00>\x00\x00\x00\x00\x00\x00\x00?\x00\x00\x00\x00\x00\x00\x00H\x00\x00\x00\x00\x00\x00\x00'
tbsS'ctx'
p14
g3
(g4
(I0
tS'b'
tRp15
(I1
(I68
tg7
I00
S'\xe9\x03\x00\x00\x00\x00\x00\x00\xea\x03\x00\x00\x00\x00\x00\x00\xeb\x03\x00\x00\x00\x00\x00\x00\xed\x03\x00\x00\x00\x00\x00\x00\xee\x03\x00\x00\x00\x00\x00\x00\xef\x03\x00\x00\x00\x00\x00\x00\xf0\x03\x00\x00\x00\x00\x00\x00\xf1\x03\x00\x00\x00\x00\x00\x00\xf2\x03\x00\x00\x00\x00\x00\x00\xf3\x03\x00\x00\x00\x00\x00\x00\xf4\x03\x00\x00\x00\x00\x00\x00\xf5\x03\x00\x00\x00\x00\x00\x00\xf6\x03\x00\x00\x00\x00\x00\x00\xf7\x03\x00\x00\x00\x00\x00\x00\xf8\x03\x00\x00\x00\x00\x00\x00\xf9\x03\x00\x00\x00\x00\x00\x00\xfa\x03\x00\x00\x00\x00\x00\x00\xfb\x03\x00\x00\x00\x00\x00\x00\xfc\x03\x00\x00\x00\x00\x00\x00\xfd\x03\x00\x00\x00\x00\x00\x00\xfe\x03\x00\x00\x00\x00\x00\x00\xff\x03\x00\x00\x00\x00\x00\x00\x00\x04\x00\x00\x00\x00\x00\x00\x01\x04\x00\x00\x00\x00\x00\x00\x02\x04\x00\x00\x00\x00\x00\x00\x03\x04\x00\x00\x00\x00\x00\x00\x04\x04\x00\x00\x00\x00\x00\x00\x05\x04\x00\x00\x00\x00\x00\x00\x06\x04\x00\x00\x00\x00\x00\x00\x07\x04\x00\x00\x00\x00\x00\x00\x08\x04\x00\x00\x00\x00\x00\x00\t\x04\x00\x00\x00\x00\x00\x00\n\x04\x00\x00\x00\x00\x00\x00\x0b\x04\x00\x00\x00\x00\x00\x00\xd1\x07\x00\x00\x00\x00\x00\x00\xd2\x07\x00\x00\x00\x00\x00\x00\xd3\x07\x00\x00\x00\x00\x00\x00\xd5\x07\x00\x00\x00\x00\x00\x00\xd6\x07\x00\x00\x00\x00\x00\x00\xd7\x07\x00\x00\x00\x00\x00\x00\xd8\x07\x00\x00\x00\x00\x00\x00\xd9\x07\x00\x00\x00\x00\x00\x00\xda\x07\x00\x00\x00\x00\x00\x00\xdb\x07\x00\x00\x00\x00\x00\x00\xdc\x07\x00\x00\x00\x00\x00\x00\xdd\x07\x00\x00\x00\x00\x00\x00\xde\x07\x00\x00\x00\x00\x00\x00\xdf\x07\x00\x00\x00\x00\x00\x00\xe0\x07\x00\x00\x00\x00\x00\x00\xe1\x07\x00\x00\x00\x00\x00\x00\xe2\x07\x00\x00\x00\x00\x00\x00\xe3\x07\x00\x00\x00\x00\x00\x00\xe4\x07\x00\x00\x00\x00\x00\x00\xe5\x07\x00\x00\x00\x00\x00\x00\xe6\x07\x00\x00\x00\x00\x00\x00\xe7\x07\x00\x00\x00\x00\x00\x00\xe8\x07\x00\x00\x00\x00\x00\x00\xe9\x07\x00\x00\x00\x00\x00\x00\xea\x07\x00\x00\x00\x00\x00\x00\xeb\x07\x00\x00\x00\x00\x00\x00\xec\x07\x00\x00\x00\x00\x00\x00\xed\x07\x00\x00\x00\x00\x00\x00\xee\x07\x00\x00\x00\x00\x00\x00\xef\x07\x00\x00\x00\x00\x00\x00\xf0\x07\x00\x00\x00\x00\x00\x00\xf1\x07\x00\x00\x00\x00\x00\x00\xf2\x07\x00\x00\x00\x00\x00\x00\xf3\x07\x00\x00\x00\x00\x00\x00'
tbsS'deep gray'
p16
(lp17
cnumpy.core.multiarray
scalar
p18
(g7
S'\n\x00\x00\x00\x00\x00\x00\x00'
tRp19
ag18
(g7
S'\x0b\x00\x00\x00\x00\x00\x00\x00'
tRp20
ag18
(g7
S'\x0c\x00\x00\x00\x00\x00\x00\x00'
tRp21
ag18
(g7
S'\r\x00\x00\x00\x00\x00\x00\x00'
tRp22
ag18
(g7
S'\x10\x00\x00\x00\x00\x00\x00\x00'
tRp23
ag18
(g7
S'\x11\x00\x00\x00\x00\x00\x00\x00'
tRp24
ag18
(g7
S'\x12\x00\x00\x00\x00\x00\x00\x00'
tRp25
ag18
(g7
S'\x1a\x00\x00\x00\x00\x00\x00\x00'
tRp26
ag18
(g7
S'\x1c\x00\x00\x00\x00\x00\x00\x00'
tRp27
ag18
(g7
S'1\x00\x00\x00\x00\x00\x00\x00'
tRp28
ag18
(g7
S'2\x00\x00\x00\x00\x00\x00\x00'
tRp29
ag18
(g7
S'3\x00\x00\x00\x00\x00\x00\x00'
tRp30
ag18
(g7
S'4\x00\x00\x00\x00\x00\x00\x00'
tRp31
ag18
(g7
S'5\x00\x00\x00\x00\x00\x00\x00'
tRp32
ag18
(g7
S'6\x00\x00\x00\x00\x00\x00\x00'
tRp33
ag18
(g7
S':\x00\x00\x00\x00\x00\x00\x00'
tRp34
ag18
(g7
S'<\x00\x00\x00\x00\x00\x00\x00'
tRp35
as.

msd-response.py

#!/usr/bin/env python
import os
import cPickle
from warnings import warn

import numpy as np
from scipy import ndimage

from dipy.reconst import shm
from dipy.reconst.dti import TensorModel
from dipy.core.geometry import vec2vec_rotmat, cart2sphere
from dipy.reconst.msd import (closest, sh_const)

_this_dir, _this_file = os.path.split(__file__)
keyfile = "{}/fs_keys.p".format(_this_dir)
fskeys = cPickle.load(open(keyfile))


def in1d(a, t):
    return np.in1d(a, t).reshape(a.shape)


def biggest(blobs):
    L, N = ndimage.label(blobs)
    sizes = np.bincount(L.ravel())
    sizes[0] = 0
    return L.reshape(blobs.shape) == sizes.argmax()


def make_wm_response_function(gtab, data, sh_order, shell_id, dirs=None,
                              no_b0=False):
    m = 0
    n = np.arange(0, sh_order + 1, 2)
    n_shells = shell_id.max() + 1

    response = np.zeros([n_shells, len(n)])
    z_unit_vector = np.array([0., 0., 1.])

    if data.ndim != 2:
        raise ValueError("Array too deep.")

    if dirs is None:
        tensor_model = TensorModel(gtab, fit_method='WLS')
        tensor_fit = tensor_model.fit(data)
        dirs = tensor_fit.evecs[..., 0]
    elif dirs.shape != (data.shape[:-1] + (3,)):
        raise ValueError("dirs has wrong shape.")

    if no_b0:
        first_dwi_shell = 0
    else:
        n_b0s = gtab.b0s_mask.sum()
        sum_b0_sig = data[:, gtab.b0s_mask].sum()
        response[0, 0] = sum_b0_sig / (sh_const * n_b0s)
        first_dwi_shell = 1

    for i_vox in range(0, data.shape[0]):
        rotmat = vec2vec_rotmat(dirs[i_vox], z_unit_vector)
        rot_gradients = np.dot(rotmat, gtab.gradients.T).T

        for i in range(first_dwi_shell, n_shells):
            shell = shell_id == i
            x, y, z = rot_gradients[shell].T
            r, theta, phi = cart2sphere(x, y, z)
            B_dwi = shm.real_sph_harm(m, n, theta[:, None], phi[:, None])
            sig_r_sh = np.linalg.lstsq(B_dwi, data[i_vox, shell])[0]
            response[i] += sig_r_sh

    response /= data.shape[0]

    return response, n


def make_isotropic_response(data, tissue_map, shell_id):
    """ Makes the isotropic part of a multi tissue response function.

    Parameters
    ----------
    data : array, shape=(..., N)
        The diffusion weighted data for fitting the response.
    tissue_map : array
        A map of voxels that belong to each tissue type.
    shell_id : array, shape=(N,)
        Which q-shell each dwi volume belongs to.

    Returns:
    --------
    response : array
        The response function.

    """
    # tissues go from 1 to n
    n_tissues = tissue_map.max()
    # shells go from 0 to n - 1
    n_shells = shell_id.max() + 1
    response = np.zeros((n_shells, n_tissues))
    for tissue in range(1, n_tissues + 1):
        for id in range(0, n_shells):
            shell = (shell_id == id)
            shell_index = shell.nonzero()[0].reshape((-1, 1))
            sub_data = data[tissue_map == tissue, shell_index]
            response[id, tissue - 1] = sub_data.mean()
            print(id, tissue - 1, response[id, tissue - 1] )
    response /= shm.real_sph_harm(0, 0, 0, 0)
    return response


def make_response(data, gtab, shells, wm_map, isotropic_tissue_map, sh_order,
                  dirs=None):
    """ Estimates multi-tissue response function from data.

    Parameters
    ----------
    data : array
        Diffusion weighted images.
    gtab : GradientTable
        The diffusion gradients used to acquire ``data``.
    shells : array
        The bvalues of each q-space shell. If the data includes non-diffusion
        weighted data, the first shell should be 0.
    wm_map : array, dtype=bool
        The voxels to use for fitting the white matter response.
    isotropic_tissue_map : array
        The voxels associated with each of the isotropic components. 0 should
        be background. Components should be represented by ints from 1 to N for
        N components.
    sh_order : int (even)
        The SH order of the white matter response.

    Returns
    -------
    response : array, shape=(n_shells, n_tissues)
        The response function as a 2d array. ``n_shells`` is the number of
        shells in the data including the b0 shell. ``n_tissues`` is the
        number of tissue types, 1 for white matter plus the number of isotropic
        tissues in the ``isotropic_tissue_map``.

    """
    shell_id = closest(shells, gtab.bvals)
    iso_response = make_isotropic_response(data, isotropic_tissue_map, shell_id)

    no_b0 = shells[0] != 0
    wm_map = np.asarray(wm_map, dtype=bool)
    wm_data = data[wm_map]
    wm_response, n = make_wm_response_function(gtab, wm_data, sh_order,
                                               shell_id, None, no_b0)
    return np.concatenate([iso_response, wm_response], axis=1)


def iso_tissue_maps_from_freesurfer(asegvol, min_vol=1000):
    cortex = in1d(asegvol, fskeys['ctx'])
    cortex = ndimage.binary_erosion(cortex, iterations=2)
    if cortex.sum() < min_vol:
        warn("Tissue volume too small, maybe something went wrong?")
    csf = in1d(asegvol, fskeys['csf'])
    csf = ndimage.binary_erosion(csf, iterations=1)
    if csf.sum() < min_vol:
        warn("Tissue volume too small, maybe something went wrong?")

    return csf + 2 * cortex


def wm_tissue_map_from_fa(fa, fa_threshold, min_vol=1000):
    wm_map = biggest(fa > fa_threshold)
    if wm_map.sum() < min_vol:
        warn("Tissue volume too small, maybe something went wrong?")
    return wm_map


def response_from_freesurfer(asegvol, data, gtab, shells, fa, sh_order):
    iso_map = iso_tissue_maps_from_freesurfer(asegvol)
    wm_map = wm_tissue_map_from_fa(fa, .7)
    response = make_response(data, gtab, shells, wm_map, iso_map, sh_order)
    return response, iso_map, wm_map