|
0
|
1 #!/usr/bin/env python
|
|
|
2 import argparse
|
|
12
|
3 import os
|
|
0
|
4 import shutil
|
|
|
5
|
|
13
|
6 from dipy.data import read_stanford_hardi
|
|
0
|
7 from dipy.reconst import peaks, shm
|
|
|
8 from dipy.tracking import utils
|
|
|
9 from dipy.tracking.eudx import EuDX
|
|
|
10 from dipy.viz import fvtk
|
|
|
11 from dipy.viz.colormap import line_colors
|
|
|
12
|
|
|
13 import matplotlib.pyplot as plt
|
|
|
14
|
|
|
15 import nibabel as nib
|
|
|
16
|
|
|
17 import numpy as np
|
|
|
18
|
|
|
19 parser = argparse.ArgumentParser()
|
|
12
|
20 parser.add_argument('--input', dest='input', help='Input dataset')
|
|
|
21 parser.add_argument('--input_extra_files_path', dest='input_extra_files_path', help='Input dataset extra files paths')
|
|
0
|
22 parser.add_argument('--output_superiorfrontal_nifti', dest='output_superiorfrontal_nifti', help='Output superiorfrontal nifti1 dataset')
|
|
12
|
23 parser.add_argument('--output_superiorfrontal_nifti_files_path', dest='output_superiorfrontal_nifti_files_path', help='Output superiorfrontal nifti1 extra files path')
|
|
0
|
24 parser.add_argument('--output_trackvis_header', dest='output_trackvis_header', help='Output superiorfrontal track visualization header dataset')
|
|
|
25
|
|
|
26 args = parser.parse_args()
|
|
|
27
|
|
12
|
28 def move_directory_files(source_dir, destination_dir, copy=False, remove_source_dir=False):
|
|
|
29 source_directory = os.path.abspath(source_dir)
|
|
|
30 destination_directory = os.path.abspath(destination_dir)
|
|
|
31 if not os.path.isdir(destination_directory):
|
|
|
32 os.makedirs(destination_directory)
|
|
|
33 for dir_entry in os.listdir(source_directory):
|
|
|
34 source_entry = os.path.join(source_directory, dir_entry)
|
|
|
35 if copy:
|
|
|
36 shutil.copy(source_entry, destination_directory)
|
|
|
37 else:
|
|
|
38 shutil.move(source_entry, destination_directory)
|
|
|
39 if remove_source_dir:
|
|
|
40 os.rmdir(source_directory)
|
|
|
41
|
|
|
42 # Get input data.
|
|
|
43 # TODO: do not hard-code 'stanford_hardi'
|
|
|
44 input_dir = 'stanford_hardi'
|
|
|
45 os.mkdir(input_dir)
|
|
|
46 for f in os.listdir(args.input_extra_files_path):
|
|
|
47 shutil.copy(os.path.join(args.input_extra_files_path, f), input_dir)
|
|
14
|
48 hardi_img, gtab = read_stanford_hardi()
|
|
0
|
49 data = hardi_img.get_data()
|
|
12
|
50 labels_file = os.path.join(input_dir, "aparc-reduced.nii.gz")
|
|
|
51 labels_img = nib.load(labels_file)
|
|
0
|
52 labels = labels_img.get_data()
|
|
|
53
|
|
|
54 white_matter = (labels == 1) | (labels == 2)
|
|
|
55 csamodel = shm.CsaOdfModel(gtab, 6)
|
|
|
56 csapeaks = peaks.peaks_from_model(model=csamodel, data=data, sphere=peaks.default_sphere, relative_peak_threshold=.8, min_separation_angle=45, mask=white_matter)
|
|
|
57 seeds = utils.seeds_from_mask(white_matter, density=2)
|
|
|
58 streamline_generator = EuDX(csapeaks.peak_values, csapeaks.peak_indices, odf_vertices=peaks.default_sphere.vertices, a_low=.05, step_sz=.5, seeds=seeds)
|
|
|
59 affine = streamline_generator.affine
|
|
|
60 streamlines = list(streamline_generator)
|
|
|
61 cc_slice = labels == 2
|
|
|
62 cc_streamlines = utils.target(streamlines, cc_slice, affine=affine)
|
|
|
63 cc_streamlines = list(cc_streamlines)
|
|
12
|
64
|
|
15
|
65 M, grouping = utils.connectivity_matrix(cc_streamlines, labels, affine=affine, return_mapping=True, mapping_as_streamlines=True)
|
|
0
|
66 lr_superiorfrontal_track = grouping[11, 54]
|
|
|
67 shape = labels.shape
|
|
|
68 dm = utils.density_map(lr_superiorfrontal_track, shape, affine=affine)
|
|
|
69 # Save density map
|
|
|
70 dm_img = nib.Nifti1Image(dm.astype("int16"), hardi_img.affine)
|
|
|
71 dm_img.to_filename("lr-superiorfrontal-dm.nii")
|
|
|
72 shutil.move('lr-superiorfrontal-dm.nii', args.output_superiorfrontal_nifti)
|
|
12
|
73 move_directory_files(input_dir, args.output_superiorfrontal_nifti_files_path)
|
|
|
74
|
|
0
|
75 # Make a trackvis header so we can save streamlines
|
|
|
76 voxel_size = labels_img.header.get_zooms()
|
|
|
77 trackvis_header = nib.trackvis.empty_header()
|
|
|
78 trackvis_header['voxel_size'] = voxel_size
|
|
|
79 trackvis_header['dim'] = shape
|
|
|
80 trackvis_header['voxel_order'] = "RAS"
|
|
|
81 # Move streamlines to "trackvis space"
|
|
|
82 trackvis_point_space = utils.affine_for_trackvis(voxel_size)
|
|
|
83 lr_sf_trk = utils.move_streamlines(lr_superiorfrontal_track, trackvis_point_space, input_space=affine)
|
|
|
84 lr_sf_trk = list(lr_sf_trk)
|
|
|
85 # Save streamlines
|
|
|
86 for_save = [(sl, None, None) for sl in lr_sf_trk]
|
|
|
87 nib.trackvis.write("lr-superiorfrontal.trk", for_save, trackvis_header)
|
|
|
88 shutil.move('lr-superiorfrontal.trk', args.output_trackvis_header)
|
