Newer
Older
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
83
84
85
86
87
88
89
90
91
#!/usr/bin/env fslpython
import os
import sys
sys.path.append('/home/paulmc/Projects/fslpy/')
from pathlib import Path
import numpy as np
from fsl.data.image import Image
from fsl.wrappers import vecreg, concatxfm, invxfm, convertwarp
from pyfeeds.evaluate import cmpVectorArrays
def main():
fsldir = Path(os.environ['FSLDIR'])
outdir = Path(sys.argv[1])
indir = Path(sys.argv[2])/'fsl_course_data'/'fdt2'
dyads = indir/'diffusion.bedpostX'/'dyads1.nii.gz'
std = fsldir/'data'/'standard'/'MNI152_T1_2mm.nii.gz'
str2std = indir/'diffusion.bedpostX'/'xfms'/'str2standard.mat'
diff2str = indir/'structural'/'xfms'/'diff2str.mat'
diff2std_warp = indir/'structural'/'xfms'/'diff2standard_warp.nii.gz'
# Create diff2std/std2diff affines
std2diff = outdir/'std2diff.mat'
diff2std = outdir/'diff2std.mat'
concatxfm(diff2str, str2std, diff2std)
invxfm(diff2std, std2diff)
# create a copy of the diff2std warp
# which does not encode the diff2std
# affine
diff2std_warp_no_premat = outdir/'diff2std_warp_no_premat.nii.gz'
convertwarp(diff2std_warp_no_premat, std,
warp1=diff2std_warp,
premat=std2diff,
relout=True)
# Transform dyads with the full diff2std
# warp (which encodes the diff->standard
# affine transform) - this will be used
# as a benchmark
dyads_std = outdir/'dyads1_std.nii.gz'
vecreg(dyads, dyads_std, std, warpfield=diff2std_warp)
# Transform dyads with the raw diff2std
# warp, and specify the diff2std affine
# as premat. The result should match
# that from the first vecreg call above.
dyads_std_with_premat = outdir/'dyads1_std_with_premat.nii.gz'
vecreg(dyads, dyads_std_with_premat, std,
warpfield=diff2std_warp_no_premat,
premat=diff2std)
# compare the result
dyads1 = Image(dyads_std).data
dyads2 = Image(dyads_std_with_premat).data
diff = cmpVectorArrays(dyads1, dyads2)
assert np.isclose(diff, 0, atol=1e-4), diff
# Create a nonsense warp field with a postmat
diff2std_warp_postmat = outdir/'diff2std_warp_postmat.nii.gz'
convertwarp(diff2std_warp_postmat, std,
warp1=diff2std_warp_no_premat,
premat=diff2std,
postmat=std2diff)
# Warp dyads with the nonsense warp
dyads_postmat1 = outdir/'dyads_postmat1.nii.gz'
vecreg(dyads, dyads_postmat1, std,
warpfield=diff2std_warp_postmat)
# And with the original warp+premat+postmat
dyads_postmat2 = outdir/'dyads_postmat2.nii.gz'
vecreg(dyads, dyads_postmat2, std,
warpfield=diff2std_warp_no_premat,
premat=diff2std, postmat=std2diff)
# And check that the results are equivalent
dyads1 = Image(dyads_postmat1).data
dyads2 = Image(dyads_postmat2).data
diff = cmpVectorArrays(dyads1, dyads2)
assert np.isclose(diff, 0, atol=1e-4), diff
if __name__ == '__main__':
main()