import platform
import dipy as dp
import matplotlib.pyplot as plt
from os.path import expanduser, join
from dipy.io import read_bvals_bvecs
from dipy.core.gradients import gradient_table
from dipy.io.image import load_nifti,save_nifti # 保存和读取nifti
print("DIPY VERSION:",dp.__version__)
print("PYTHON VERSION:",platform.python_version())

DIPY VERSION: 1.5.0
PYTHON VERSION: 3.8.5

dMRI(diffusion MRI)经常使用三种类型的文件，包含扩散像权重的NIfTI数据、包含b-values的文本文件(.bval)、包含b-vectors的文本文件(.bvec).

# 下载一个dMRI数据集，在~/.dipy/sherbrooke_3shell目录内。
from dipy.data import fetch_sherbrooke_3shell
fetch_sherbrooke_3shell()

({'HARDI193.nii.gz': ('https://digital.lib.washington.edu/researchworks/bitstream/handle/1773/38475/HARDI193.nii.gz',
   '0b735e8f16695a37bfbd66aab136eb66'),
  'HARDI193.bval': ('https://digital.lib.washington.edu/researchworks/bitstream/handle/1773/38475/HARDI193.bval',
   'e9b9bb56252503ea49d31fb30a0ac637'),
  'HARDI193.bvec': ('https://digital.lib.washington.edu/researchworks/bitstream/handle/1773/38475/HARDI193.bvec',
   '0c83f7e8b917cd677ad58a078658ebb7')},
 '/Users/reallo/.dipy/sherbrooke_3shell')

# 获得home directory
home_path = expanduser('~')
home_path

'/Users/reallo'

dname = join(home_path, '.dipy', 'sherbrooke_3shell')
dname

'/Users/reallo/.dipy/sherbrooke_3shell'

fdwi = join(dname, 'HARDI193.nii.gz')

print(fdwi)

fbval = join(dname, 'HARDI193.bval')

print(fbval)

fbvec = join(dname, 'HARDI193.bvec')

print(fbvec)

/Users/reallo/.dipy/sherbrooke_3shell/HARDI193.nii.gz
/Users/reallo/.dipy/sherbrooke_3shell/HARDI193.bval
/Users/reallo/.dipy/sherbrooke_3shell/HARDI193.bvec

加载dMRI数据集

data,affine,img = load_nifti(fdwi,return_img=True)

print("data数据的大小：",data.shape)
print("dimensions of each voxel:",img.header.get_zooms()[:3])

data数据的大小： (128, 128, 60, 193)
dimensions of each voxel: (2.0, 2.0, 2.0)

该数据为multi volume文件，包含了193个volume。 该图使用MRIcron软件查看

affine

array([[   2.        ,    0.        ,   -0.        , -128.02410889],
       [   0.        ,    2.        ,   -0.        ,  -97.02124023],
       [  -0.        ,   -0.        ,    2.        ,  -73.61566162],
       [   0.        ,    0.        ,    0.        ,    1.        ]])

type(data)

numpy.ndarray

data[:,:,:,0]

array([[[12, 24, 15, ...,  0,  0,  0],
        [27, 15, 19, ..., 15, 15, 16],
        [12, 24, 13, ..., 17, 13, 20],
        ...,
        [ 9, 12, 20, ..., 18, 23, 20],
        [20, 17, 11, ..., 22, 15, 12],
        [29, 18, 12, ..., 16, 11, 14]],

       [[11, 19, 15, ...,  0,  0,  0],
        [14, 12, 16, ..., 31, 18, 22],
        [11, 21, 17, ..., 15, 20, 15],
        ...,
        [26, 26, 33, ..., 17, 17, 24],
        [15,  8, 27, ..., 10, 16, 15],
        [25,  6, 22, ..., 19, 17, 22]],

       [[18,  8, 15, ...,  0,  0,  0],
        [21, 22, 29, ..., 20, 16, 22],
        [18,  7, 16, ..., 20, 15, 22],
        ...,
        [27, 34, 16, ..., 21, 14, 24],
        [25, 15, 38, ..., 20, 18, 24],
        [16, 27, 29, ..., 19, 14, 25]],

       ...,

       [[22,  7, 16, ...,  0,  0,  0],
        [20, 24, 27, ..., 28, 19, 27],
        [22,  4, 16, ..., 24, 15, 16],
        ...,
        [13, 16, 26, ..., 17, 37,  9],
        [16, 26, 24, ..., 22, 31, 23],
        [20, 21, 16, ..., 19, 26, 20]],

       [[25, 18, 31, ...,  0,  0,  0],
        [25, 30, 26, ..., 16, 18, 25],
        [25, 17, 35, ...,  8, 27, 19],
        ...,
        [21, 15, 20, ..., 22, 25, 24],
        [24, 23, 25, ..., 21, 29, 32],
        [31, 26, 27, ..., 19, 27, 24]],

       [[ 0,  0,  0, ...,  0,  0,  0],
        [ 0,  0,  0, ...,  0,  0,  0],
        [ 0,  0,  0, ...,  0,  0,  0],
        ...,
        [ 0,  0,  0, ...,  0,  0,  0],
        [ 0,  0,  0, ...,  0,  0,  0],
        [ 0,  0,  0, ...,  0,  0,  0]]], dtype=uint16)

img

<nibabel.nifti1.Nifti1Image at 0x102ca26d0>

使用matplotlib库画图

axial_middle = data.shape[2]//2 
axial_middle 

30

plt.figure()
plt.subplot(1,2,1).set_axis_off() # 设置子图(左)
plt.imshow(data[:,:,axial_middle,0].T,cmap='gray',origin='lower')
plt.title('volume=0,axial slice')

plt.subplot(1,2,2).set_axis_off() # 设置子图(右)
plt.title('volume=10,diffusion weighting')
plt.imshow(data[:,:,axial_middle,10].T,cmap='gray',origin='lower')

<matplotlib.image.AxesImage at 0x11e3b66d0>

plt.close()

读取bval和bvec文件

bvals,bvecs = read_bvals_bvecs(fbval,fbvec)

print("bvals shape:",bvals.shape)
bvals

bvals shape: (193,)

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., 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., 2000., 2000., 2000., 2000., 2000., 2000., 2000.,
       2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000.,
       2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000.,
       2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000.,
       2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000.,
       2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000.,
       2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000.,
       2000., 2000., 2000., 3500., 3500., 3500., 3500., 3500., 3500.,
       3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500.,
       3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500.,
       3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500.,
       3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500.,
       3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500.,
       3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500.,
       3500., 3500., 3500., 3500.])

print("bvecs.shape:",bvecs.shape)
bvecs[0:10,:]

bvecs.shape: (193, 3)

array([[ 0.        ,  0.        ,  0.        ],
       [ 0.999979  , -0.00504001, -0.00402795],
       [ 0.        ,  0.999992  , -0.00398794],
       [-0.0257055 ,  0.653861  , -0.756178  ],
       [ 0.589518  , -0.769236  , -0.246462  ],
       [-0.235785  , -0.529095  , -0.815147  ],
       [-0.893578  , -0.263559  , -0.363394  ],
       [ 0.79784   ,  0.133726  , -0.587851  ],
       [ 0.232937  ,  0.931884  , -0.278087  ],
       [ 0.93672   ,  0.144139  , -0.31903   ]])

#### 使用GradientTable来记录采集时的特定参数，包括b-values,b-vectors,timing等

gtab = gradient_table(bvals,bvecs)

gtab.info

B-values shape (193,)
         min 0.000000
         max 3500.000000
B-vectors shape (193, 3)
          min -0.964050
          max 0.999992

# 查看b-values 
gtab.bvals

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., 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., 2000., 2000., 2000., 2000., 2000., 2000., 2000.,
       2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000.,
       2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000.,
       2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000.,
       2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000.,
       2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000.,
       2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000., 2000.,
       2000., 2000., 2000., 3500., 3500., 3500., 3500., 3500., 3500.,
       3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500.,
       3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500.,
       3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500.,
       3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500.,
       3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500.,
       3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500., 3500.,
       3500., 3500., 3500., 3500.])

dir(gtab)

['__class__',
 '__delattr__',
 '__dict__',
 '__dir__',
 '__doc__',
 '__eq__',
 '__format__',
 '__ge__',
 '__getattribute__',
 '__gt__',
 '__hash__',
 '__init__',
 '__init_subclass__',
 '__le__',
 '__lt__',
 '__module__',
 '__ne__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__sizeof__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 'b0_threshold',
 'b0s_mask',
 'big_delta',
 'btens',
 'bvals',
 'bvecs',
 'gradient_strength',
 'gradients',
 'info',
 'qvals',
 'small_delta',
 'tau']

gtab.b0_threshold

50

## 前10个b-vectors
gtab.bvecs[:10,:]

array([[ 0.        ,  0.        ,  0.        ],
       [ 0.999979  , -0.00504001, -0.00402795],
       [ 0.        ,  0.999992  , -0.00398794],
       [-0.0257055 ,  0.653861  , -0.756178  ],
       [ 0.589518  , -0.769236  , -0.246462  ],
       [-0.235785  , -0.529095  , -0.815147  ],
       [-0.893578  , -0.263559  , -0.363394  ],
       [ 0.79784   ,  0.133726  , -0.587851  ],
       [ 0.232937  ,  0.931884  , -0.278087  ],
       [ 0.93672   ,  0.144139  , -0.31903   ]])

gtab可以告诉我们S0 volumes(b-values=0)

gtab.b0s_mask

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, 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, 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, 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, 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, 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, False, False, False, False, False, False, False, False,
       False, False, False, False])

S0s = data[:,:,:,gtab.b0s_mask]
S0s.shape

(128, 128, 60, 1)

保存为一个新的NIfTI文件

save_nifti("HARDI193_S0.nii.gz",S0s,affine) # 保存到当前文件夹