Mercurial > repos > greg > create_streamlines

comparison create_streamlines.py @ 12:883327c42970 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Thu, 30 Nov 2017 10:11:00 -0500
parents 1f1fdfe9ac4d
children ab72c2d00e5a
comparison
equal deleted inserted replaced
11:7b98048e1904 12:883327c42970
1 #!/usr/bin/env python 1 #!/usr/bin/env python
2 import argparse 2 import argparse
3 import os
3 import shutil 4 import shutil
4 5
5 from dipy.data import fetch_stanford_t1, read_stanford_labels, read_stanford_t1
6 from dipy.reconst import peaks, shm 6 from dipy.reconst import peaks, shm
7 from dipy.tracking import utils 7 from dipy.tracking import utils
8 from dipy.tracking.eudx import EuDX 8 from dipy.tracking.eudx import EuDX
9 from dipy.viz import fvtk 9 from dipy.viz import fvtk
10 from dipy.viz.colormap import line_colors 10 from dipy.viz.colormap import line_colors
14 import nibabel as nib 14 import nibabel as nib
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('--input', dest='input', help='Input dataset')
20 parser.add_argument('--input_extra_files_path', dest='input_extra_files_path', help='Input dataset extra files paths')
20 parser.add_argument('--output_superiorfrontal_nifti', dest='output_superiorfrontal_nifti', help='Output superiorfrontal nifti1 dataset') 21 parser.add_argument('--output_superiorfrontal_nifti', dest='output_superiorfrontal_nifti', help='Output superiorfrontal nifti1 dataset')
22 parser.add_argument('--output_superiorfrontal_nifti_files_path', dest='output_superiorfrontal_nifti_files_path', help='Output superiorfrontal nifti1 extra files path')
21 parser.add_argument('--output_trackvis_header', dest='output_trackvis_header', help='Output superiorfrontal track visualization header dataset') 23 parser.add_argument('--output_trackvis_header', dest='output_trackvis_header', help='Output superiorfrontal track visualization header dataset')
22 24
23 args = parser.parse_args() 25 args = parser.parse_args()
24 26
25 hardi_img, gtab, labels_img = read_stanford_labels() 27 def move_directory_files(source_dir, destination_dir, copy=False, remove_source_dir=False):
28 source_directory = os.path.abspath(source_dir)
29 destination_directory = os.path.abspath(destination_dir)
30 if not os.path.isdir(destination_directory):
31 os.makedirs(destination_directory)
32 for dir_entry in os.listdir(source_directory):
33 source_entry = os.path.join(source_directory, dir_entry)
34 if copy:
35 shutil.copy(source_entry, destination_directory)
36 else:
37 shutil.move(source_entry, destination_directory)
38 if remove_source_dir:
39 os.rmdir(source_directory)
40
41 # Get input data.
42 # TODO: do not hard-code 'stanford_hardi'
43 input_dir = 'stanford_hardi'
44 os.mkdir(input_dir)
45 for f in os.listdir(args.input_extra_files_path):
46 shutil.copy(os.path.join(args.input_extra_files_path, f), input_dir)
47 hard_img, gtab = read_stanford_hardi()
26 data = hardi_img.get_data() 48 data = hardi_img.get_data()
49 labels_file = os.path.join(input_dir, "aparc-reduced.nii.gz")
50 labels_img = nib.load(labels_file)
27 labels = labels_img.get_data() 51 labels = labels_img.get_data()
28 52
29 # For possible later use: if args.drmi_dataset == 'stanford_t1':
30 #fetch_stanford_t1()
31 #t1 = read_stanford_t1()
32 #t1_data = t1.get_data()
33 white_matter = (labels == 1) | (labels == 2) 53 white_matter = (labels == 1) | (labels == 2)
34 csamodel = shm.CsaOdfModel(gtab, 6) 54 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) 55 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) 56 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) 57 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 58 affine = streamline_generator.affine
39 streamlines = list(streamline_generator) 59 streamlines = list(streamline_generator)
40 cc_slice = labels == 2 60 cc_slice = labels == 2
41 cc_streamlines = utils.target(streamlines, cc_slice, affine=affine) 61 cc_streamlines = utils.target(streamlines, cc_slice, affine=affine)
42 cc_streamlines = list(cc_streamlines) 62 cc_streamlines = list(cc_streamlines)
43 #other_streamlines = utils.target(streamlines, cc_slice, affine=affine, include=False) 63
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)
66 M, grouping = utils.connectivity_matrix(cc_streamlines, labels, affine=affine, return_mapping=True, mapping_as_streamlines=True)
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)
72 lr_superiorfrontal_track = grouping[11, 54] 64 lr_superiorfrontal_track = grouping[11, 54]
73 shape = labels.shape 65 shape = labels.shape
74 dm = utils.density_map(lr_superiorfrontal_track, shape, affine=affine) 66 dm = utils.density_map(lr_superiorfrontal_track, shape, affine=affine)
75 # Save density map 67 # Save density map
76 dm_img = nib.Nifti1Image(dm.astype("int16"), hardi_img.affine) 68 dm_img = nib.Nifti1Image(dm.astype("int16"), hardi_img.affine)
77 dm_img.to_filename("lr-superiorfrontal-dm.nii") 69 dm_img.to_filename("lr-superiorfrontal-dm.nii")
78 shutil.move('lr-superiorfrontal-dm.nii', args.output_superiorfrontal_nifti) 70 shutil.move('lr-superiorfrontal-dm.nii', args.output_superiorfrontal_nifti)
71 move_directory_files(input_dir, args.output_superiorfrontal_nifti_files_path)
72
79 # Make a trackvis header so we can save streamlines 73 # Make a trackvis header so we can save streamlines
80 voxel_size = labels_img.header.get_zooms() 74 voxel_size = labels_img.header.get_zooms()
81 trackvis_header = nib.trackvis.empty_header() 75 trackvis_header = nib.trackvis.empty_header()
82 trackvis_header['voxel_size'] = voxel_size 76 trackvis_header['voxel_size'] = voxel_size
83 trackvis_header['dim'] = shape 77 trackvis_header['dim'] = shape