-
Notifications
You must be signed in to change notification settings - Fork 1
/
write_fg_to_trk_shift.m
148 lines (127 loc) · 4.82 KB
/
write_fg_to_trk_shift.m
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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
function write_fg_to_trk_shift( fg, nii_file, trk_file)
% Write fg data into Trackvis format
%
% write_fg_to_trk(fg, filename)
%
% write_fg_to_trk write fiber group structure fg to trackvis file.
%
% Input:
% fg: is a mrDiffusion fiber group structure
% nii_file: the nifti image for reference space of fibers
% trk_file: Name of trk file to output
% Recommended the use of '.trk' extension of file
% Output:
%
%
% For details about header fields and fileformat see:
% http://www.trackvis.org/docs/?subsect=fileformat
%
% Example;
%
% write_fg_to_trk(fg, 'dwi.nii.gz', 'csd_life.trk');
%
% Load the header of nifti reference file
hdr_nii = niftiRead(nii_file);
% Dummy empty header
hdr_trk.id_string = 'TRACK';
% The dimensions of the reference nifti file
hdr_trk.dim = hdr_nii.dim(1:3);
% The size of the voxels
hdr_trk.voxel_size = hdr_nii.pixdim(1:3);
% Apparently not implemented by Trackvis
% hdr_trk.origin = [hdr_nii.qoffset_x, hdr_nii.qoffset_y, hdr_nii.qoffset_z];
hdr_trk.origin = [0.0, 0.0, 0.0];
% Optional features not managed by the current implementation
hdr_trk.n_scalars = 0;
hdr_trk.scalar_name = ['','','','','','','','','',''];
hdr_trk.n_properties = 0;
hdr_trk.property_name = ['','','','','','','','','',''];
% The affine transformation
scale_dim = diag(1 ./ hdr_trk.voxel_size);
affine = hdr_nii.qto_xyz;
% correct displacement of Vistasoft
affine(1,4) = affine(1,4) + affine(1,1);
affine(2,4) = affine(2,4) + affine(2,2);
affine(3,4) = affine(3,4) + affine(3,3);
hdr_trk.vox_to_ras = transpose(affine);
affine(1:3,1:3) = scale_dim .* affine(1:3,1:3);
affine_trackvis = affine;
affine_trackvis(1,4) = affine_trackvis(1,4) + hdr_trk.voxel_size(1)/2;
affine_trackvis(2:3,4) = affine_trackvis(2:3,4) - hdr_trk.voxel_size(1)/2;
nii2trk = inv(affine_trackvis);
hdr_trk.reserved = '';
% Original orientation of Trackvis toolbox
hdr_trk.voxel_order = 'LAS';
% Orientation of the reference image
hdr_trk.pad2 = '';
hdr_trk.image_orientation_patient = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
hdr_trk.pad1 = '';
hdr_trk.invert_x = '';
hdr_trk.invert_y = '';
hdr_trk.invert_z = '';
hdr_trk.swap_xy = '';
hdr_trk.swap_yz = '';
hdr_trk.swap_zx = '';
% Total number of fibers/streamlines
hdr_trk.n_count = size(fg.fibers,1);
% Version of the current header
hdr_trk.version = 2;
% The size of the header as number of byte
hdr_trk.hdr_size = 1000;
% Begin writing the trk file
fid = fopen(trk_file , 'w');
% Writing the header of trk file
% (some features are not implemented)
skip = 0;
wb = fwrite(fid,hdr_trk.id_string,'char');
skip = 6-wb*1;
wb = fwrite(fid,hdr_trk.dim,'3*int16',skip);
skip = 6-wb*2 + skip * (1 - sign(wb));
wb = fwrite(fid,hdr_trk.voxel_size,'3*float',skip);
skip = 12-wb*4 + skip * (1 - sign(wb));
wb = fwrite(fid,hdr_trk.origin,'3*float',skip);
skip = 12-wb*4 + skip * (1 - sign(wb));
wb = fwrite(fid,hdr_trk.n_scalars,'int16',skip);
skip = 2-wb*2 + skip * (1 - sign(wb));
wb = fwrite(fid,hdr_trk.scalar_name,'200*char',skip);
skip = 200-wb*1 + skip * (1 - sign(wb));
wb = fwrite(fid,hdr_trk.n_properties,'int16',skip);
skip = 2-wb*2 + skip * (1 - sign(wb));
wb = fwrite(fid,hdr_trk.property_name,'200*char',skip);
skip = 200-wb*1 + skip * (1 - sign(wb));
wb = fwrite(fid,hdr_trk.vox_to_ras,'16*float',skip);
skip = 64-wb*4 + skip * (1 - sign(wb));
wb = fwrite(fid,hdr_trk.reserved,'444*char',skip);
skip = 444-wb*1 + skip * (1 - sign(wb));
wb = fwrite(fid,hdr_trk.voxel_order,'4*char',skip);
skip = 4-wb*1 + skip * (1 - sign(wb));
wb = fwrite(fid,hdr_trk.pad2,'4*char',skip);
skip = 4-wb*1 + skip * (1 - sign(wb));
wb = fwrite(fid,hdr_trk.image_orientation_patient,'6*float',skip);
skip = 24-wb*4 + skip * (1 - sign(wb));
wb = fwrite(fid, hdr_trk.pad1,'2*char', skip);
skip = 2-wb*1 + skip * (1 - sign(wb));
wb = fwrite(fid, hdr_trk.invert_x,'1*uchar',skip);
skip = 1-wb*1 + skip * (1 - sign(wb));
wb = fwrite(fid, hdr_trk.invert_y,'1*uchar', skip);
skip = 1-wb*1 + skip * (1 - sign(wb));
wb = fwrite(fid, hdr_trk.invert_z,'1*uchar', skip);skip = 1-wb*1 + skip * (1 - sign(wb));
wb = fwrite(fid, hdr_trk.swap_xy,'1*uchar', skip);
skip = 1-wb*1 + skip * (1 - sign(wb));
wb = fwrite(fid, hdr_trk.swap_yz,'1*uchar', skip);
skip = 1-wb*1 + skip * (1 - sign(wb));
wb = fwrite(fid, hdr_trk.swap_zx,'1*uchar', skip);
skip = 1-wb*1 + skip * (1 - sign(wb));
wb = fwrite(fid, hdr_trk.n_count,'1*int32', skip);
skip = 4-wb*4 + skip * (1 - sign(wb));
wb = fwrite(fid, hdr_trk.version,'1*int32', skip);
skip = 4-wb*4 + skip * (1 - sign(wb));
wb = fwrite(fid, hdr_trk.hdr_size,'1*int32', skip);
% Writing the fibers/streamlines points to trk file
% (currently not supported the optional scalar and track properties)
for t=1:size(fg.fibers,1)
fwrite(fid, size(fg.fibers{t},2),'int32');
fpoints = mrAnatXformCoords(nii2trk, transpose(fg.fibers{t}));
fwrite(fid,reshape(transpose(fpoints),1,numel(fpoints)),'float');
end
fclose(fid)