-
Notifications
You must be signed in to change notification settings - Fork 0
/
tractograms_slr.py
82 lines (63 loc) · 2.98 KB
/
tractograms_slr.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
""" SLR (Streamline Linear Registration) of two tractograms.
See Garyfallidis et. al, 'Robust and efficient linear registration
of white-matter fascicles in the space of streamlines',
Neuroimage, 117:124-140, 2015.
"""
import os
import sys
import argparse
import nibabel as nib
import numpy as np
import ntpath
from os.path import isfile
from shutil import copyfile
from nibabel.streamlines import load
from dipy.segment.clustering import QuickBundles
from dipy.align.streamlinear import StreamlineLinearRegistration
from dipy.tracking.streamline import set_number_of_points
def tractograms_slr(moving_tractogram, static_tractogram):
subjID = ntpath.basename(static_tractogram)[0:6]
exID = ntpath.basename(moving_tractogram)[0:6]
aff_dir = '/N/dc2/projects/lifebid/giulia/data/HCP3-IU-Giulia/derivatives/slr_transformations'
affine_path = '%s/affine_m%s_s%s.npy' %(aff_dir, exID, subjID)
affine_fname = './affine_m%s_s%s.npy' %(exID, subjID)
if isfile(affine_path):
print("Affine already computed. Retrieving past results.")
copyfile(affine_path, affine_fname)
else:
print("Loading tractograms...")
moving_tractogram = nib.streamlines.load(moving_tractogram)
moving_tractogram = moving_tractogram.streamlines
static_tractogram = nib.streamlines.load(static_tractogram)
static_tractogram = static_tractogram.streamlines
print("Set parameters as in Garyfallidis et al. 2015.")
threshold_length = 40.0 # 50mm / 1.25
qb_threshold = 16.0 # 20mm / 1.25
nb_res_points = 20
print("Performing QuickBundles of static tractogram and resampling...")
st = np.array([s for s in static_tractogram if len(s) > threshold_length], dtype=np.object)
qb = QuickBundles(threshold=qb_threshold)
st_clusters = [cluster.centroid for cluster in qb.cluster(st)]
st_clusters = set_number_of_points(st_clusters, nb_res_points)
print("Performing QuickBundles of moving tractogram and resampling...")
mt = np.array([s for s in moving_tractogram if len(s) > threshold_length], dtype=np.object)
qb = QuickBundles(threshold=qb_threshold)
mt_clusters = [cluster.centroid for cluster in qb.cluster(mt)]
mt_clusters = set_number_of_points(mt_clusters, nb_res_points)
print("Performing Linear Registration...")
srr = StreamlineLinearRegistration()
srm = srr.optimize(static=st_clusters, moving=mt_clusters)
print("Affine transformation matrix with Streamline Linear Registration:")
affine = srm.matrix
print('%s' %affine)
np.save('affine_m%s_s%s.npy' %(exID, subjID), affine)
print("Affine for example %s and target %s saved." %(exID, subjID))
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument('-moving', nargs='?', const=1, default='',
help='The moving tractogram filename')
parser.add_argument('-static', nargs='?', const=1, default='',
help='The static tractogram filename')
args = parser.parse_args()
tractograms_slr(args.moving, args.static)
sys.exit()