|
0
|
1 #!/usr/bin/env python
|
|
|
2 import argparse
|
|
9
|
3 import os
|
|
0
|
4 import shutil
|
|
|
5
|
|
|
6 from dipy.data import fetch_sherbrooke_3shell
|
|
|
7 from dipy.data import fetch_stanford_hardi
|
|
|
8 from dipy.data import get_sphere
|
|
|
9 from dipy.data import read_sherbrooke_3shell
|
|
|
10 from dipy.data import read_stanford_hardi
|
|
|
11 from dipy.reconst import dti
|
|
|
12 from dipy.reconst.dti import color_fa
|
|
|
13 from dipy.reconst.dti import fractional_anisotropy
|
|
|
14 from dipy.segment.mask import median_otsu
|
|
|
15 from dipy.viz import fvtk
|
|
|
16 from matplotlib import pyplot
|
|
|
17
|
|
|
18 import nibabel
|
|
|
19 import numpy
|
|
|
20
|
|
|
21 # http://nipy.org/dipy/examples_built/reconst_dti.html#example-reconst-dti
|
|
|
22 parser = argparse.ArgumentParser()
|
|
9
|
23 parser.add_argument('--input', dest='input', help='Input dataset')
|
|
|
24 parser.add_argument('--input_extra_files_path', dest='input_extra_files_path', help='Input dataset extra files path')
|
|
0
|
25 parser.add_argument('--output_nifti1_fa', dest='output_nifti1_fa', help='Output fractional anisotropy Nifti1 dataset')
|
|
9
|
26 parser.add_argument('--output_nifti1_fa_files_path', dest='output_nifti1_fa_files_path', help='Output fractional anisotropy Nifti1 extra files path')
|
|
3
|
27 parser.add_argument('--output_nifti1_evecs', dest='output_nifti1_evecs', help='Output eigen vectors Nifti1 dataset')
|
|
9
|
28 parser.add_argument('--output_nifti1_evecs_files_path', dest='output_nifti1_evecs_files_path', help='Output eigen vectors Nifti1 extra files path')
|
|
0
|
29 parser.add_argument('--output_nifti1_md', dest='output_nifti1_md', help='Output mean diffusivity Nifti1 dataset')
|
|
9
|
30 parser.add_argument('--output_nifti1_md_files_path', dest='output_nifti1_md_files_path', help='Output mean diffusivity Nifti1 extra files path')
|
|
0
|
31 parser.add_argument('--output_nifti1_rgb', dest='output_nifti1_rgb', help='Output RGB-map Nifti1 dataset')
|
|
9
|
32 parser.add_argument('--output_nifti1_rgb_files_path', dest='output_nifti1_rgb_files_path', help='Output RGB-map Nifti1 extra files path')
|
|
0
|
33 parser.add_argument('--output_png_ellipsoids', dest='output_png_ellipsoids', help='Output ellipsoids PNG dataset')
|
|
|
34 parser.add_argument('--output_png_odfs', dest='output_png_odfs', help='Output orientation distribution functions PNG dataset')
|
|
|
35 parser.add_argument('--output_png_middle_axial_slice', dest='output_png_middle_axial_slice', help='Output middle axial slice PNG dataset')
|
|
|
36
|
|
|
37 args = parser.parse_args()
|
|
|
38
|
|
9
|
39 def move_directory_files(source_dir, destination_dir, copy=False, remove_source_dir=False):
|
|
|
40 source_directory = os.path.abspath(source_dir)
|
|
|
41 destination_directory = os.path.abspath(destination_dir)
|
|
|
42 if not os.path.isdir(destination_directory):
|
|
|
43 os.makedirs(destination_directory)
|
|
|
44 for dir_entry in os.listdir(source_directory):
|
|
|
45 source_entry = os.path.join(source_directory, dir_entry)
|
|
|
46 if copy:
|
|
|
47 shutil.copy(source_entry, destination_directory)
|
|
|
48 else:
|
|
|
49 shutil.move(source_entry, destination_directory)
|
|
|
50 if remove_source_dir:
|
|
|
51 os.rmdir(source_directory)
|
|
|
52
|
|
0
|
53 # Get input data.
|
|
9
|
54 # TODO: do not hard-code 'stanford_hardi'
|
|
|
55 input_dir = 'stanford_hardi'
|
|
|
56 os.mkdir(input_dir)
|
|
0
|
57 input_dir = args.drmi_dataset
|
|
9
|
58 for f in os.list_dir(args.input_extra_files_path):
|
|
|
59 shutil.copy(f, input_dir)
|
|
|
60 fetch_stanford_hardi()
|
|
|
61 img, gtab = read_stanford_hardi()
|
|
0
|
62
|
|
|
63 data = img.get_data()
|
|
|
64 maskdata, mask = median_otsu(data, 3, 1, True, vol_idx=range(10, 50), dilate=2)
|
|
|
65
|
|
|
66 axial_middle = data.shape[2] // 2
|
|
|
67 pyplot.subplot(1, 2, 1).set_axis_off()
|
|
|
68 pyplot.imshow(data[:, :, axial_middle, 0].T, cmap='gray', origin='lower')
|
|
|
69 pyplot.subplot(1, 2, 2).set_axis_off()
|
|
|
70 pyplot.imshow(data[:, :, axial_middle, 10].T, cmap='gray', origin='lower')
|
|
|
71 pyplot.savefig('middle_axial.png', bbox_inches='tight')
|
|
|
72 shutil.move('middle_axial.png', args.output_png_middle_axial_slice)
|
|
|
73
|
|
|
74 tenmodel = dti.TensorModel(gtab)
|
|
|
75 tenfit = tenmodel.fit(maskdata)
|
|
|
76
|
|
|
77 fa = fractional_anisotropy(tenfit.evals)
|
|
|
78 fa[numpy.isnan(fa)] = 0
|
|
|
79 fa_img = nibabel.Nifti1Image(fa.astype(numpy.float32), img.affine)
|
|
|
80 nibabel.save(fa_img, 'output_fa.nii')
|
|
|
81 shutil.move('output_fa.nii', args.output_nifti1_fa)
|
|
9
|
82 move_directory_files(input_dir, args.output_nifti1_fa_files_path)
|
|
0
|
83
|
|
|
84 evecs_img = nibabel.Nifti1Image(tenfit.evecs.astype(numpy.float32), img.affine)
|
|
|
85 nibabel.save(evecs_img, 'output_evecs.nii')
|
|
|
86 shutil.move('output_evecs.nii', args.output_nifti1_evecs)
|
|
9
|
87 move_directory_files(input_dir, args.output_nifti1_evecs_files_path)
|
|
0
|
88
|
|
|
89 md1 = dti.mean_diffusivity(tenfit.evals)
|
|
|
90 nibabel.save(nibabel.Nifti1Image(md1.astype(numpy.float32), img.affine), 'output_md.nii')
|
|
|
91 shutil.move('output_md.nii', args.output_nifti1_md)
|
|
9
|
92 move_directory_files(input_dir, args.output_nifti1_md_files_path)
|
|
0
|
93
|
|
|
94 fa = numpy.clip(fa, 0, 1)
|
|
|
95 rgb = color_fa(fa, tenfit.evecs)
|
|
|
96 nibabel.save(nibabel.Nifti1Image(numpy.array(255 * rgb, 'uint8'), img.affine), 'output_rgb.nii')
|
|
|
97 shutil.move('output_rgb.nii', args.output_nifti1_rgb)
|
|
9
|
98 move_directory_files(input_dir, args.output_nifti1_rgb_files_path)
|
|
0
|
99
|
|
|
100 sphere = get_sphere('symmetric724')
|
|
|
101 ren = fvtk.ren()
|
|
|
102
|
|
|
103 evals = tenfit.evals[13:43, 44:74, 28:29]
|
|
|
104 evecs = tenfit.evecs[13:43, 44:74, 28:29]
|
|
|
105 cfa = rgb[13:43, 44:74, 28:29]
|
|
|
106 cfa /= cfa.max()
|
|
|
107 fvtk.add(ren, fvtk.tensor(evals, evecs, cfa, sphere))
|
|
|
108 fvtk.record(ren, n_frames=1, out_path='tensor_ellipsoids.png', size=(600, 600))
|
|
|
109 shutil.move('tensor_ellipsoids.png', args.output_png_ellipsoids)
|
|
|
110
|
|
|
111 fvtk.clear(ren)
|
|
|
112
|
|
|
113 tensor_odfs = tenmodel.fit(data[20:50, 55:85, 38:39]).odf(sphere)
|
|
|
114 fvtk.add(ren, fvtk.sphere_funcs(tensor_odfs, sphere, colormap=None))
|
|
|
115 fvtk.record(ren, n_frames=1, out_path='tensor_odfs.png', size=(600, 600))
|
|
|
116 shutil.move('tensor_odfs.png', args.output_png_odfs)
|
