-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathread_ctf_hc.m
185 lines (169 loc) · 7.82 KB
/
read_ctf_hc.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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
function [hc] = read_ctf_hc(filename);
% READ_CTF_HC reads the MEG headcoil marker positions from an ascii file
% and computes the coordinate transformation required to get from from
% dewar to head-coordinates
%
% the definition of head coordinates is according to CTF standard:
% - the origin is exactly between LPA and RPA
% - the positive x-axis goes throught NAS
% - the positive y-axis goes (approximately) through LPA
% - the positive z-axis goes up, orthogonal to the x- and y-axes
%
% hc = read_ctf_hc(filename)
%
% returns a structure with the following fields
% hc.dewar.nas marker positions relative to dewar
% hc.dewar.lpa
% hc.dewar.rpa
% hc.head.nas marker positions relative to head (measured)
% hc.head.lpa
% hc.head.rpa
% hc.standard.nas marker positions relative to head (expected)
% hc.standard.lpa
% hc.standard.rpa
% and
% hc.affine parameter for affine transformation (1x12)
% hc.homogenous homogenous transformation matrix (4x4, see warp3d)
% hc.translation translation vector (1x3)
% hc.rotation rotation matrix (3x3)
%
% Gradiometer positions can be transformed into head coordinates using the
% homogeneous transformation matrix, or using the affine parameters and
% the warp3d function from the WARPING toolbox
%
% WARNING: all returned positions are in mm, and the transformation matrices
% also should be applied on channel positions that are in mm.
% Copyright (C) 2002, Robert Oostenveld
%
global fb
hc.standard.nas = [0 0 0];
hc.standard.lpa = [0 0 0];
hc.standard.rpa = [0 0 0];
hc.dewar.nas = [0 0 0];
hc.dewar.lpa = [0 0 0];
hc.dewar.rpa = [0 0 0];
hc.head.nas = [0 0 0];
hc.head.lpa = [0 0 0];
hc.head.rpa = [0 0 0];
fid = fopen(filename, 'r');
if fid==-1
error(sprintf('could not open file %s', filename));
end
fseek(fid, 0, 'bof');
line = [];
if fb
fprintf('reading standard coil positions with respect to the dewar\n');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% standard coil positions with respect to the dewar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while ~strcmp(line, 'standard nasion coil position relative to dewar (cm):')
line = fgetl(fid);
if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.nas(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.nas(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.nas(3) = str2num(r(2:end));
% NOTE THAT THERE IS AN TYPING ERROR IN SOME CTF FILES WHICH I HAVE TO REPRODUCE HERE (staNdard)
while ~(strcmp(line, 'stadard left ear coil position relative to dewar (cm):') | ...
strcmp(line, 'standard left ear coil position relative to dewar (cm):'))
line = fgetl(fid);
if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.lpa(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.lpa(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.lpa(3) = str2num(r(2:end));
while ~strcmp(line, 'standard right ear coil position relative to dewar (cm):')
line = fgetl(fid);
if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.rpa(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.rpa(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.rpa(3) = str2num(r(2:end));
if fb
fprintf('reading measured coil positions with respect to the dewar\n');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% measured coil positions with respect to the dewar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while ~strcmp(line, 'measured nasion coil position relative to dewar (cm):')
line = fgetl(fid);
if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.nas(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.nas(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.nas(3) = str2num(r(2:end));
while ~strcmp(line, 'measured left ear coil position relative to dewar (cm):')
line = fgetl(fid);
if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.lpa(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.lpa(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.lpa(3) = str2num(r(2:end));
while ~strcmp(line, 'measured right ear coil position relative to dewar (cm):')
line = fgetl(fid);
if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.rpa(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.rpa(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.rpa(3) = str2num(r(2:end));
if fb
fprintf('reading measured coil positions with respect to the head\n');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% measured coil positions with respect to the head
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while ~strcmp(line, 'measured nasion coil position relative to head (cm):')
line = fgetl(fid);
if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.nas(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.nas(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.nas(3) = str2num(r(2:end));
while ~strcmp(line, 'measured left ear coil position relative to head (cm):')
line = fgetl(fid);
if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.lpa(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.lpa(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.lpa(3) = str2num(r(2:end));
while ~strcmp(line, 'measured right ear coil position relative to head (cm):')
line = fgetl(fid);
if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.rpa(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.rpa(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.rpa(3) = str2num(r(2:end));
fclose(fid);
if fb
fprintf('computing transformation matrix\n');
end
% change the marker positions from cm to mm
hc.standard.nas = 10*hc.standard.nas;
hc.standard.lpa = 10*hc.standard.lpa;
hc.standard.rpa = 10*hc.standard.rpa;
hc.dewar.nas = 10*hc.dewar.nas;
hc.dewar.lpa = 10*hc.dewar.lpa;
hc.dewar.rpa = 10*hc.dewar.rpa;
hc.head.nas = 10*hc.head.nas;
hc.head.lpa = 10*hc.head.lpa;
hc.head.rpa = 10*hc.head.rpa;
% compute the direction of the head coordinate axes in dewar coordinates
d_x = hc.dewar.nas - (hc.dewar.lpa + hc.dewar.rpa)/2;
d_z = cross(d_x, hc.dewar.lpa - hc.dewar.rpa);
d_y = cross(d_z, d_x);
d_x = d_x / norm(d_x);
d_y = d_y / norm(d_y);
d_z = d_z / norm(d_z);
% compute the translation and rotation which are neccessary to transform
% any given location from the dewar to the head coordinate system
hc.translation = -(hc.dewar.lpa + hc.dewar.rpa)/2;
hc.rotation = inv([d_x' d_y' d_z']);
% construct the homogenous coordinate transformation matrix
% from the subsequent translation and rotation
hc.homogenous = eye(4,4);
hc.homogenous(1:3,4) = hc.rotation * hc.translation';
hc.homogenous(1:3,1:3) = hc.rotation;
% construct a vector with the 12 parameters for an affine warp
hc.affine = hc.homogenous';
hc.affine = hc.affine(1:12);