|
0
|
1 #!/usr/bin/env python
|
|
|
2 import argparse
|
|
|
3 import shutil
|
|
|
4
|
|
|
5 from dipy.data import fetch_stanford_t1, read_stanford_labels, read_stanford_t1
|
|
|
6 from dipy.reconst import peaks, shm
|
|
|
7 from dipy.tracking import utils
|
|
|
8 from dipy.tracking.eudx import EuDX
|
|
|
9 from dipy.viz import fvtk
|
|
|
10 from dipy.viz.colormap import line_colors
|
|
|
11
|
|
|
12 import matplotlib.pyplot as plt
|
|
|
13
|
|
|
14 import nibabel as nib
|
|
|
15
|
|
|
16 import numpy as np
|
|
|
17
|
|
|
18 parser = argparse.ArgumentParser()
|
|
|
19 parser.add_argument('--drmi_dataset', dest='drmi_dataset', help='Input dataset')
|
|
|
20 parser.add_argument('--output_superiorfrontal_nifti', dest='output_superiorfrontal_nifti', help='Output superiorfrontal nifti1 dataset')
|
|
|
21 parser.add_argument('--output_trackvis_header', dest='output_trackvis_header', help='Output superiorfrontal track visualization header dataset')
|
|
|
22
|
|
|
23 args = parser.parse_args()
|
|
|
24
|
|
|
25 hardi_img, gtab, labels_img = read_stanford_labels()
|
|
|
26 data = hardi_img.get_data()
|
|
|
27 labels = labels_img.get_data()
|
|
|
28
|
|
3
|
29 # For possible later use: if args.drmi_dataset == 'stanford_t1':
|
|
8
|
30 #fetch_stanford_t1()
|
|
|
31 #t1 = read_stanford_t1()
|
|
|
32 #t1_data = t1.get_data()
|
|
0
|
33 white_matter = (labels == 1) | (labels == 2)
|
|
|
34 csamodel = shm.CsaOdfModel(gtab, 6)
|
|
|
35 csapeaks = peaks.peaks_from_model(model=csamodel, data=data, sphere=peaks.default_sphere, relative_peak_threshold=.8, min_separation_angle=45, mask=white_matter)
|
|
|
36 seeds = utils.seeds_from_mask(white_matter, density=2)
|
|
|
37 streamline_generator = EuDX(csapeaks.peak_values, csapeaks.peak_indices, odf_vertices=peaks.default_sphere.vertices, a_low=.05, step_sz=.5, seeds=seeds)
|
|
|
38 affine = streamline_generator.affine
|
|
|
39 streamlines = list(streamline_generator)
|
|
|
40 cc_slice = labels == 2
|
|
|
41 cc_streamlines = utils.target(streamlines, cc_slice, affine=affine)
|
|
|
42 cc_streamlines = list(cc_streamlines)
|
|
8
|
43 #other_streamlines = utils.target(streamlines, cc_slice, affine=affine, include=False)
|
|
|
44 #other_streamlines = list(other_streamlines)
|
|
|
45 #assert len(other_streamlines) + len(cc_streamlines) == len(streamlines)
|
|
|
46 ## Make display objects
|
|
|
47 #color = line_colors(cc_streamlines)
|
|
|
48 #cc_streamlines_actor = fvtk.line(cc_streamlines, line_colors(cc_streamlines))
|
|
|
49 #cc_ROI_actor = fvtk.contour(cc_slice, levels=[1], colors=[(1., 1., 0.)], opacities=[1.])
|
|
|
50 #vol_actor = fvtk.slicer(t1_data)
|
|
|
51 #vol_actor.display(40, None, None)
|
|
|
52 #vol_actor2 = vol_actor.copy()
|
|
|
53 #vol_actor2.display(None, None, 35)
|
|
|
54 ## Add display objects to canvas
|
|
|
55 #r = fvtk.ren()
|
|
|
56 #fvtk.add(r, vol_actor)
|
|
|
57 #fvtk.add(r, vol_actor2)
|
|
|
58 #fvtk.add(r, cc_streamlines_actor)
|
|
|
59 #fvtk.add(r, cc_ROI_actor)
|
|
|
60 ## Save figures
|
|
|
61 #fvtk.record(r, n_frames=1, out_path="corpuscallosum_axial.png", size=(800, 800))
|
|
|
62 #shutil.move("corpuscallosum_axial.png", args.output_corpuscallosum_axial)
|
|
|
63 #fvtk.camera(r, [-1, 0, 0], [0, 0, 0], viewup=[0, 0, 1])
|
|
|
64 #fvtk.record(r, n_frames=1, out_path="corpuscallosum_sagittal.png", size=(800, 800))
|
|
|
65 #shutil.move("corpuscallosum_sagittal.png", args.output_corpuscallosum_sagittal)
|
|
0
|
66 M, grouping = utils.connectivity_matrix(cc_streamlines, labels, affine=affine, return_mapping=True, mapping_as_streamlines=True)
|
|
8
|
67 #M[:3, :] = 0
|
|
|
68 #M[:, :3] = 0
|
|
|
69 #plt.imshow(np.log1p(M), interpolation='nearest')
|
|
|
70 #plt.savefig("connectivity.png")
|
|
|
71 #shutil.move("connectivity.png", args.output_connectivity)
|
|
0
|
72 lr_superiorfrontal_track = grouping[11, 54]
|
|
|
73 shape = labels.shape
|
|
|
74 dm = utils.density_map(lr_superiorfrontal_track, shape, affine=affine)
|
|
|
75 # Save density map
|
|
|
76 dm_img = nib.Nifti1Image(dm.astype("int16"), hardi_img.affine)
|
|
|
77 dm_img.to_filename("lr-superiorfrontal-dm.nii")
|
|
|
78 shutil.move('lr-superiorfrontal-dm.nii', args.output_superiorfrontal_nifti)
|
|
|
79 # Make a trackvis header so we can save streamlines
|
|
|
80 voxel_size = labels_img.header.get_zooms()
|
|
|
81 trackvis_header = nib.trackvis.empty_header()
|
|
|
82 trackvis_header['voxel_size'] = voxel_size
|
|
|
83 trackvis_header['dim'] = shape
|
|
|
84 trackvis_header['voxel_order'] = "RAS"
|
|
|
85 # Move streamlines to "trackvis space"
|
|
|
86 trackvis_point_space = utils.affine_for_trackvis(voxel_size)
|
|
|
87 lr_sf_trk = utils.move_streamlines(lr_superiorfrontal_track, trackvis_point_space, input_space=affine)
|
|
|
88 lr_sf_trk = list(lr_sf_trk)
|
|
|
89 # Save streamlines
|
|
|
90 for_save = [(sl, None, None) for sl in lr_sf_trk]
|
|
|
91 nib.trackvis.write("lr-superiorfrontal.trk", for_save, trackvis_header)
|
|
|
92 shutil.move('lr-superiorfrontal.trk', args.output_trackvis_header)
|
