Python dti.TensorModel类代码示例

def test_eudx_bad_seed():
    """Test passing a bad seed to eudx"""
    fimg, fbvals, fbvecs = get_data('small_101D')

    img = ni.load(fimg)
    affine = img.get_affine()
    data = img.get_data()
    gtab = gradient_table(fbvals, fbvecs)
    tensor_model = TensorModel(gtab)
    ten = tensor_model.fit(data)
    ind = quantize_evecs(ten.evecs)

    seed = [1000000., 1000000., 1000000.]
    eu = EuDX(a=ten.fa, ind=ind, seeds=[seed], a_low=.2)
    try:
        track = list(eu)
    except ValueError as ve:        
        if ve.args[0] == 'Seed outside boundaries':
            print(ve)
   
    print(data.shape)
    seed = [1., 5., 8.]
    eu = EuDX(a=ten.fa, ind=ind, seeds=[seed], a_low=.2)    
    track = list(eu)
    
    seed = [-1., 1000000., 1000000.]
    eu = EuDX(a=ten.fa, ind=ind, seeds=[seed], a_low=.2)
    try:
        track = list(eu)
    except ValueError as ve:
        if ve.args[0] == 'Seed outside boundaries':
            print(ve)

def test_eudx_bad_seed():
    """Test passing a bad seed to eudx"""
    fimg, fbvals, fbvecs = get_data('small_101D')

    img = ni.load(fimg)
    affine = img.affine
    data = img.get_data()
    gtab = gradient_table(fbvals, fbvecs)
    tensor_model = TensorModel(gtab)
    ten = tensor_model.fit(data)
    ind = quantize_evecs(ten.evecs)

    sphere = get_sphere('symmetric724')
    seed = [1000000., 1000000., 1000000.]
    eu = EuDX(a=ten.fa, ind=ind, seeds=[seed],
              odf_vertices=sphere.vertices, a_low=.2)
    assert_raises(ValueError, list, eu)

    print(data.shape)
    seed = [1., 5., 8.]
    eu = EuDX(a=ten.fa, ind=ind, seeds=[seed],
              odf_vertices=sphere.vertices, a_low=.2)
    track = list(eu)

    seed = [-1., 1000000., 1000000.]
    eu = EuDX(a=ten.fa, ind=ind, seeds=[seed],
              odf_vertices=sphere.vertices, a_low=.2)
    assert_raises(ValueError, list, eu)

def test_eudx_further():
    """ Cause we love testin.. ;-)
    """

    fimg,fbvals,fbvecs=get_data('small_101D')

    img=ni.load(fimg)
    affine=img.get_affine()
    data=img.get_data()
    gtab = gradient_table(fbvals, fbvecs)
    tensor_model = TensorModel(gtab)
    ten = tensor_model.fit(data)
    x,y,z=data.shape[:3]
    seeds=np.zeros((10**4,3))
    for i in range(10**4):
        rx=(x-1)*np.random.rand()
        ry=(y-1)*np.random.rand()
        rz=(z-1)*np.random.rand()            
        seeds[i]=np.ascontiguousarray(np.array([rx,ry,rz]),dtype=np.float64)
    
    ind = quantize_evecs(ten.evecs)
    eu=EuDX(a=ten.fa, ind=ind, seeds=seeds, a_low=.2)
    T=[e for e in eu]
    
    #check that there are no negative elements
    for t in T:
        assert_equal(np.sum(t.ravel()<0),0)

def test_response_from_mask():
    fdata, fbvals, fbvecs = get_data('small_64D')
    bvals = np.load(fbvals)
    bvecs = np.load(fbvecs)
    data = nib.load(fdata).get_data()

    gtab = gradient_table(bvals, bvecs)
    ten = TensorModel(gtab)
    tenfit = ten.fit(data)
    FA = fractional_anisotropy(tenfit.evals)
    FA[np.isnan(FA)] = 0
    radius = 3

    for fa_thr in np.arange(0, 1, 0.1):
        response_auto, ratio_auto, nvoxels = auto_response(gtab,
                                                           data,
                                                           roi_center=None,
                                                           roi_radius=radius,
                                                           fa_thr=fa_thr,
                                                           return_number_of_voxels=True)

        ci, cj, ck = np.array(data.shape[:3]) / 2
        mask = np.zeros(data.shape[:3])
        mask[ci - radius: ci + radius,
             cj - radius: cj + radius,
             ck - radius: ck + radius] = 1

        mask[FA <= fa_thr] = 0
        response_mask, ratio_mask = response_from_mask(gtab, data, mask)

        assert_equal(int(np.sum(mask)), nvoxels)
        assert_array_almost_equal(response_mask[0], response_auto[0])
        assert_almost_equal(response_mask[1], response_auto[1])
        assert_almost_equal(ratio_mask, ratio_auto)

def compute_tensor_model(dir_src, dir_out, verbose=False):

    fbval = pjoin(dir_src, 'bvals_' + par_b_tag)
    fbvec = pjoin(dir_src, 'bvecs_' + par_b_tag)
    fdwi =  pjoin(dir_src, 'data_' + par_b_tag + '_' + par_dim_tag + '.nii.gz')
    fmask = pjoin(dir_src, 'nodif_brain_mask_' + par_dim_tag + '.nii.gz')

    bvals, bvecs = read_bvals_bvecs(fbval, fbvec)
    gtab = gradient_table(bvals, bvecs, b0_threshold=par_b0_threshold)
    data, affine = load_nifti(fdwi, verbose)
    mask, _ = load_nifti(fmask, verbose)

    ten_model = TensorModel(gtab)
    ten_fit = ten_model.fit(data, mask)

    FA = ten_fit.fa
    MD = ten_fit.md
    EV = ten_fit.evecs.astype(np.float32)

    fa_name = 'data_' + par_b_tag + '_' + par_dim_tag + '_FA.nii.gz'
    save_nifti(pjoin(dir_out, fa_name), FA, affine)
    md_name = 'data_' + par_b_tag + '_' + par_dim_tag + '_MD.nii.gz'
    save_nifti(pjoin(dir_out, md_name), MD, affine)
    ev_name = 'data_' + par_b_tag + '_' + par_dim_tag + '_EV.nii.gz'
    save_nifti(pjoin(dir_out, ev_name), EV, affine)

def test_phantom():
    N = 50

    vol = orbital_phantom(gtab,
                          func=f,
                          t=np.linspace(0, 2 * np.pi, N),
                          datashape=(10, 10, 10, len(bvals)),
                          origin=(5, 5, 5),
                          scale=(3, 3, 3),
                          angles=np.linspace(0, 2 * np.pi, 16),
                          radii=np.linspace(0.2, 2, 6),
                          S0=100)

    m = TensorModel(gtab)
    t = m.fit(vol)
    FA = t.fa
    # print vol
    FA[np.isnan(FA)] = 0
    # 686 -> expected FA given diffusivities of [1500, 400, 400]
    l1, l2, l3 = 1500e-6, 400e-6, 400e-6
    expected_fa = (np.sqrt(0.5) *
                   np.sqrt((l1 - l2)**2 + (l2-l3)**2 + (l3-l1)**2) /
                   np.sqrt(l1**2 + l2**2 + l3**2))

    assert_array_almost_equal(FA.max(), expected_fa, decimal=2)

def test_masked_array_with_tensor():
    data = np.ones((2, 4, 56))
    mask = np.array([[True, False, False, True],
                     [True, False, True, False]])

    bvec, bval = read_bvec_file(get_data('55dir_grad.bvec'))
    gtab = grad.gradient_table_from_bvals_bvecs(bval, bvec.T)

    tensor_model = TensorModel(gtab)
    tensor = tensor_model.fit(data, mask=mask)
    assert_equal(tensor.shape, (2, 4))
    assert_equal(tensor.fa.shape, (2, 4))
    assert_equal(tensor.evals.shape, (2, 4, 3))
    assert_equal(tensor.evecs.shape, (2, 4, 3, 3))

    tensor = tensor[0]
    assert_equal(tensor.shape, (4,))
    assert_equal(tensor.fa.shape, (4,))
    assert_equal(tensor.evals.shape, (4, 3))
    assert_equal(tensor.evecs.shape, (4, 3, 3))

    tensor = tensor[0]
    assert_equal(tensor.shape, tuple())
    assert_equal(tensor.fa.shape, tuple())
    assert_equal(tensor.evals.shape, (3,))
    assert_equal(tensor.evecs.shape, (3, 3))
    assert_equal(type(tensor.model_params), np.ndarray)

def DIPY_nii2streamlines(imgfile, maskfile, bvals, bvecs, output_prefix):
    import numpy as np
    import nibabel as nib
    import os

    from dipy.reconst.dti import TensorModel

    print "nii2streamlines"

    img = nib.load(imgfile)
    bvals = np.genfromtxt(bvals)
    bvecs = np.genfromtxt(bvecs)
    if bvecs.shape[1] != 3:
        bvecs = bvecs.T
    print bvecs.shape

    from nipype.utils.filemanip import split_filename
    _, prefix, _  = split_filename(imgfile)
    from dipy.data import gradient_table
    gtab = gradient_table(bvals, bvecs)
    data = img.get_data()
    affine = img.get_affine()
    zooms = img.get_header().get_zooms()[:3]
    new_zooms = (2., 2., 2.)
    data2, affine2 = data, affine
    mask = nib.load(maskfile).get_data().astype(np.bool)
    tenmodel = TensorModel(gtab)
    tenfit = tenmodel.fit(data2, mask)

    from dipy.reconst.dti import fractional_anisotropy
    FA = fractional_anisotropy(tenfit.evals)
    FA[np.isnan(FA)] = 0
    fa_img = nib.Nifti1Image(FA, img.get_affine())
    nib.save(fa_img, experiment_dir + '/' + ('%s_tensor_fa.nii.gz' % prefix))
    evecs = tenfit.evecs
    evec_img = nib.Nifti1Image(evecs, img.get_affine())
    nib.save(evec_img, experiment_dir + '/' + ('%s_tensor_evec.nii.gz' % prefix))

    from dipy.data import get_sphere
    sphere = get_sphere('symmetric724')
    from dipy.reconst.dti import quantize_evecs

    peak_indices = quantize_evecs(tenfit.evecs, sphere.vertices)

    from dipy.tracking.eudx import EuDX
    eu = EuDX(FA, peak_indices, odf_vertices = sphere.vertices, a_low=0.2, seeds=10**6, ang_thr=35)
    tensor_streamlines = [streamline for streamline in eu]
    hdr = nib.trackvis.empty_header()
    hdr['voxel_size'] = new_zooms
    hdr['voxel_order'] = 'LPS'
    hdr['dim'] = data2.shape[:3]

    import dipy.tracking.metrics as dmetrics
    tensor_streamlines = ((sl, None, None) for sl in tensor_streamlines if dmetrics.length(sl) > 15)
    ten_sl_fname = experiment_dir + '/' + ('%s_streamline.trk' % prefix)
    nib.trackvis.write(ten_sl_fname, tensor_streamlines, hdr, points_space='voxel')
    return ten_sl_fname

def test_eudx_further():
    """ Cause we love testin.. ;-)
    """

    fimg, fbvals, fbvecs = get_data('small_101D')

    img = ni.load(fimg)
    affine = img.affine
    data = img.get_data()
    gtab = gradient_table(fbvals, fbvecs)
    tensor_model = TensorModel(gtab)
    ten = tensor_model.fit(data)
    x, y, z = data.shape[:3]
    seeds = np.zeros((10**4, 3))
    for i in range(10**4):
        rx = (x-1)*np.random.rand()
        ry = (y-1)*np.random.rand()
        rz = (z-1)*np.random.rand()
        seeds[i] = np.ascontiguousarray(np.array([rx, ry, rz]),
                                        dtype=np.float64)

    sphere = get_sphere('symmetric724')

    ind = quantize_evecs(ten.evecs)
    eu = EuDX(a=ten.fa, ind=ind, seeds=seeds,
              odf_vertices=sphere.vertices, a_low=.2)
    T = [e for e in eu]

    # check that there are no negative elements
    for t in T:
        assert_equal(np.sum(t.ravel() < 0), 0)

    # Test eudx with affine
    def random_affine(seeds):
        affine = np.eye(4)
        affine[:3, :] = np.random.random((3, 4))
        seeds = np.dot(seeds, affine[:3, :3].T)
        seeds += affine[:3, 3]
        return affine, seeds

    # Make two random affines and move seeds
    affine1, seeds1 = random_affine(seeds)
    affine2, seeds2 = random_affine(seeds)

    # Make tracks using different affines
    eu1 = EuDX(a=ten.fa, ind=ind, odf_vertices=sphere.vertices,
               seeds=seeds1, a_low=.2, affine=affine1)
    eu2 = EuDX(a=ten.fa, ind=ind, odf_vertices=sphere.vertices,
               seeds=seeds2, a_low=.2, affine=affine2)

    # Move from eu2 affine2 to affine1
    eu2_to_eu1 = utils.move_streamlines(eu2, output_space=affine1,
                                        input_space=affine2)
    # Check that the tracks are the same
    for sl1, sl2 in zip(eu1, eu2_to_eu1):
        assert_array_almost_equal(sl1, sl2)

def test_single_tensor():
    evals = np.array([1.4, .35, .35]) * 10**(-3)
    evecs = np.eye(3)
    S = SingleTensor(gtab, 100, evals, evecs, snr=None)
    assert_array_almost_equal(S[gtab.b0s_mask], 100)
    assert_(np.mean(S[~gtab.b0s_mask]) < 100)

    from dipy.reconst.dti import TensorModel
    m = TensorModel(gtab)
    t = m.fit(S)

    assert_array_almost_equal(t.fa, 0.707, decimal=3)

def FA_RGB(data, gtab):
    """
    Input : data, gtab taken from the load_data.py script.
    Return : FA and RGB as two nd numpy array
    """

    tenmodel = TensorModel(gtab)
    tenfit = tenmodel.fit(data)
    FA = fractional_anisotropy(tenfit.evals)
    FA[np.isnan(FA)] = 0
    FA = np.clip(FA, 0, 1)
    RGB = color_fa(FA, tenfit.evecs)
    return FA, RGB

def estimate_response(gtab, data, affine, mask, fa_thr=0.7):
    tenmodel = TensorModel(gtab)
    tenfit = tenmodel.fit(data, mask)
    FA = fractional_anisotropy(tenfit.evals)
    FA[np.isnan(FA)] = 0
    mask[FA <= 0.1] = 0
    mask[FA > 1.] = 0
    indices = np.where(FA > fa_thr)
    lambdas = tenfit.evals[indices][:, :2]
    S0s = data[indices][:, 0]
    S0 = np.mean(S0s)
    l01 = np.mean(lambdas, axis=0)
    evals = np.array([l01[0], l01[1], l01[1]])
    ratio = evals[1] / evals[0]
    print 'Response evals' , evals, ' ratio: ', ratio, '\tMean S0', S0
    return (evals, S0), ratio

def prepare(training, category, snr, denoised, odeconv, tv, method):

    data, affine, gtab = get_specific_data(training,
                                           category,
                                           snr,
                                           denoised)

    prefix = create_file_prefix(training,
                                category,
                                snr,
                                denoised,
                                odeconv,
                                tv,
                                method)

    if training:
        mask = nib.load('wm_mask_hardi_01.nii.gz').get_data()
    else:
        #mask = np.ones(data.shape[:-1])
        mask = nib.load('test_hardi_30_den=1_fa_0025_dilate2_mask.nii.gz').get_data()

    tenmodel = TensorModel(gtab)

    tenfit = tenmodel.fit(data, mask)

    FA = fractional_anisotropy(tenfit.evals)
    FA[np.isnan(FA)] = 0

    mask[FA <= 0.1] = 0
    mask[FA > 1.] = 0

    indices = np.where(FA > 0.7)
    lambdas = tenfit.evals[indices][:, :2]
    S0s = data[indices][:, 0]
    S0 = np.mean(S0s)

    if S0 == 0:
        print 'S0 equals to 0 switching to 1'
        S0 = 1

    l01 = np.mean(lambdas, axis=0)

    evals = np.array([l01[0], l01[1], l01[1]])

    print evals, S0

    return data, affine, gtab, mask, evals, S0, prefix

def reconstruction(dwi,bval_file,bvec_file,mask=None,type='dti',b0=0.,order=4):
	""" Uses Dipy to reconstruct an fODF for each voxel.
	    
	    Parameters
 		----------
 	    dwi: numpy array (mandatory)
	    	Holds the diffusion weighted image in a 4D-array (see nibabel).
	    bval_file: string (mandatory)
	    	Path to the b-value file (FSL format).
	    bvec_file: string (mandatory)
	    	Path to the b-vectors file (FSL format).
	    mask:  numpy array
	    	Holds the mask in a 3D array (see nibabel).
	    type: string \in {'dti','csd','csa'} (default = 'dti')
	    	The type of the ODF reconstruction.
		b0: float (default = 0)
			Threshold to use for defining b0 images.
	    order: int (default = 4)
	    	Order to use for constrained spherical deconvolution (csd) or constant solid angle (csa).
	    	
	    Returns
		-----------
		model_fit: Dipy Object (depends on the type)
			Represents the fitted model for each voxel.
	"""	
	
	#b-values and b-vectors
	bvals, bvecs = read_bvals_bvecs(bval_file,bvec_file)
	gtab = gradient_table(bvals, bvecs, b0_threshold=b0)
	
	#reconstruction
	if type == 'dti':
		model = TensorModel(gtab,fit_method='WLS')
	elif type == 'csd':
		response, ratio = auto_response(gtab, dwi, roi_radius=10, fa_thr=0.7)
		model = ConstrainedSphericalDeconvModel(gtab, response, sh_order=order)
	elif type == 'csa':
		model = CsaOdfModel(gtab, order)

	if mask is not None:
		model_fit = model.fit(dwi,mask=mask)
	else:
		model_fit = model.fit(dwi)
	
	return model_fit

def response_from_mask(gtab, data, mask):
    """ Estimate the response function from a given mask.

    Parameters
    ----------
    gtab : GradientTable
    data : ndarray
        Diffusion data
    mask : ndarray
        Mask to use for the estimation of the response function. For example a
        mask of the white matter voxels with FA values higher than 0.7
        (see [1]_).

    Returns
    -------
    response : tuple, (2,)
        (`evals`, `S0`)
    ratio : float
        The ratio between smallest versus largest eigenvalue of the response.

    Notes
    -----
    See csdeconv.auto_response() or csdeconv.recursive_response() if you don't
    have a computed mask for the response function estimation.

    References
    ----------
    .. [1] Tournier, J.D., et al. NeuroImage 2004. Direct estimation of the
    fiber orientation density function from diffusion-weighted MRI
    data using spherical deconvolution
    """

    ten = TensorModel(gtab)
    indices = np.where(mask > 0)

    if indices[0].size == 0:
        msg = "No voxel in mask with value > 0 were found."
        warnings.warn(msg, UserWarning)
        return (np.nan, np.nan), np.nan

    tenfit = ten.fit(data[indices])
    lambdas = tenfit.evals[:, :2]
    S0s = data[indices][:, np.nonzero(gtab.b0s_mask)[0]]

    return _get_response(S0s, lambdas)