Mercurial > repos > greg > create_streamlines

comparison create_streamlines.py @ 8:1f1fdfe9ac4d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Wed, 29 Nov 2017 13:43:15 -0500
parents 606a712b6322
children 883327c42970
comparison
equal deleted inserted replaced
7:2d06191d7820 8:1f1fdfe9ac4d
15 15
16 import numpy as np 16 import numpy as np
17 17
18 parser = argparse.ArgumentParser() 18 parser = argparse.ArgumentParser()
19 parser.add_argument('--drmi_dataset', dest='drmi_dataset', help='Input dataset') 19 parser.add_argument('--drmi_dataset', dest='drmi_dataset', help='Input dataset')
20 parser.add_argument('--output_corpuscallosum_axial', dest='output_corpuscallosum_axial', help='Output corpuscallosum axial dataset')
21 parser.add_argument('--output_corpuscallosum_sagittal', dest='output_corpuscallosum_sagittal', help='Output corpuscallosum sagittal dataset')
22 parser.add_argument('--output_connectivity', dest='output_connectivity', help='Output connectivity dataset')
23 parser.add_argument('--output_superiorfrontal_nifti', dest='output_superiorfrontal_nifti', help='Output superiorfrontal nifti1 dataset') 20 parser.add_argument('--output_superiorfrontal_nifti', dest='output_superiorfrontal_nifti', help='Output superiorfrontal nifti1 dataset')
24 parser.add_argument('--output_trackvis_header', dest='output_trackvis_header', help='Output superiorfrontal track visualization header dataset') 21 parser.add_argument('--output_trackvis_header', dest='output_trackvis_header', help='Output superiorfrontal track visualization header dataset')
25 22
26 args = parser.parse_args() 23 args = parser.parse_args()
27 24
28 hardi_img, gtab, labels_img = read_stanford_labels() 25 hardi_img, gtab, labels_img = read_stanford_labels()
29 data = hardi_img.get_data() 26 data = hardi_img.get_data()
30 labels = labels_img.get_data() 27 labels = labels_img.get_data()
31 28
32 # For possible later use: if args.drmi_dataset == 'stanford_t1': 29 # For possible later use: if args.drmi_dataset == 'stanford_t1':
33 fetch_stanford_t1() 30 #fetch_stanford_t1()
34 t1 = read_stanford_t1() 31 #t1 = read_stanford_t1()
35 t1_data = t1.get_data() 32 #t1_data = t1.get_data()
36 white_matter = (labels == 1) | (labels == 2) 33 white_matter = (labels == 1) | (labels == 2)
37 csamodel = shm.CsaOdfModel(gtab, 6) 34 csamodel = shm.CsaOdfModel(gtab, 6)
38 csapeaks = peaks.peaks_from_model(model=csamodel, data=data, sphere=peaks.default_sphere, relative_peak_threshold=.8, min_separation_angle=45, mask=white_matter) 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)
39 seeds = utils.seeds_from_mask(white_matter, density=2) 36 seeds = utils.seeds_from_mask(white_matter, density=2)
40 streamline_generator = EuDX(csapeaks.peak_values, csapeaks.peak_indices, odf_vertices=peaks.default_sphere.vertices, a_low=.05, step_sz=.5, seeds=seeds) 37 streamline_generator = EuDX(csapeaks.peak_values, csapeaks.peak_indices, odf_vertices=peaks.default_sphere.vertices, a_low=.05, step_sz=.5, seeds=seeds)
41 affine = streamline_generator.affine 38 affine = streamline_generator.affine
42 streamlines = list(streamline_generator) 39 streamlines = list(streamline_generator)
43 cc_slice = labels == 2 40 cc_slice = labels == 2
44 cc_streamlines = utils.target(streamlines, cc_slice, affine=affine) 41 cc_streamlines = utils.target(streamlines, cc_slice, affine=affine)
45 cc_streamlines = list(cc_streamlines) 42 cc_streamlines = list(cc_streamlines)
46 other_streamlines = utils.target(streamlines, cc_slice, affine=affine, include=False) 43 #other_streamlines = utils.target(streamlines, cc_slice, affine=affine, include=False)
47 other_streamlines = list(other_streamlines) 44 #other_streamlines = list(other_streamlines)
48 assert len(other_streamlines) + len(cc_streamlines) == len(streamlines) 45 #assert len(other_streamlines) + len(cc_streamlines) == len(streamlines)
49 # Make display objects 46 ## Make display objects
50 color = line_colors(cc_streamlines) 47 #color = line_colors(cc_streamlines)
51 cc_streamlines_actor = fvtk.line(cc_streamlines, line_colors(cc_streamlines)) 48 #cc_streamlines_actor = fvtk.line(cc_streamlines, line_colors(cc_streamlines))
52 cc_ROI_actor = fvtk.contour(cc_slice, levels=[1], colors=[(1., 1., 0.)], opacities=[1.]) 49 #cc_ROI_actor = fvtk.contour(cc_slice, levels=[1], colors=[(1., 1., 0.)], opacities=[1.])
53 vol_actor = fvtk.slicer(t1_data) 50 #vol_actor = fvtk.slicer(t1_data)
54 vol_actor.display(40, None, None) 51 #vol_actor.display(40, None, None)
55 vol_actor2 = vol_actor.copy() 52 #vol_actor2 = vol_actor.copy()
56 vol_actor2.display(None, None, 35) 53 #vol_actor2.display(None, None, 35)
57 # Add display objects to canvas 54 ## Add display objects to canvas
58 r = fvtk.ren() 55 #r = fvtk.ren()
59 fvtk.add(r, vol_actor) 56 #fvtk.add(r, vol_actor)
60 fvtk.add(r, vol_actor2) 57 #fvtk.add(r, vol_actor2)
61 fvtk.add(r, cc_streamlines_actor) 58 #fvtk.add(r, cc_streamlines_actor)
62 fvtk.add(r, cc_ROI_actor) 59 #fvtk.add(r, cc_ROI_actor)
63 # Save figures 60 ## Save figures
64 fvtk.record(r, n_frames=1, out_path="corpuscallosum_axial.png", size=(800, 800)) 61 #fvtk.record(r, n_frames=1, out_path="corpuscallosum_axial.png", size=(800, 800))
65 shutil.move("corpuscallosum_axial.png", args.output_corpuscallosum_axial) 62 #shutil.move("corpuscallosum_axial.png", args.output_corpuscallosum_axial)
66 fvtk.camera(r, [-1, 0, 0], [0, 0, 0], viewup=[0, 0, 1]) 63 #fvtk.camera(r, [-1, 0, 0], [0, 0, 0], viewup=[0, 0, 1])
67 fvtk.record(r, n_frames=1, out_path="corpuscallosum_sagittal.png", size=(800, 800)) 64 #fvtk.record(r, n_frames=1, out_path="corpuscallosum_sagittal.png", size=(800, 800))
68 shutil.move("corpuscallosum_sagittal.png", args.output_corpuscallosum_sagittal) 65 #shutil.move("corpuscallosum_sagittal.png", args.output_corpuscallosum_sagittal)
69 M, grouping = utils.connectivity_matrix(cc_streamlines, labels, affine=affine, return_mapping=True, mapping_as_streamlines=True) 66 M, grouping = utils.connectivity_matrix(cc_streamlines, labels, affine=affine, return_mapping=True, mapping_as_streamlines=True)
70 M[:3, :] = 0 67 #M[:3, :] = 0
71 M[:, :3] = 0 68 #M[:, :3] = 0
72 plt.imshow(np.log1p(M), interpolation='nearest') 69 #plt.imshow(np.log1p(M), interpolation='nearest')
73 plt.savefig("connectivity.png") 70 #plt.savefig("connectivity.png")
74 shutil.move("connectivity.png", args.output_connectivity) 71 #shutil.move("connectivity.png", args.output_connectivity)
75 lr_superiorfrontal_track = grouping[11, 54] 72 lr_superiorfrontal_track = grouping[11, 54]
76 shape = labels.shape 73 shape = labels.shape
77 dm = utils.density_map(lr_superiorfrontal_track, shape, affine=affine) 74 dm = utils.density_map(lr_superiorfrontal_track, shape, affine=affine)
78 # Save density map 75 # Save density map
79 dm_img = nib.Nifti1Image(dm.astype("int16"), hardi_img.affine) 76 dm_img = nib.Nifti1Image(dm.astype("int16"), hardi_img.affine)