-
Notifications
You must be signed in to change notification settings - Fork 0
/
ctf_head2mri.m
187 lines (135 loc) · 6.16 KB
/
ctf_head2mri.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
185
186
187
function MRIVertices = ctf_head2mri(HeadVertices,mri)
% ctf_head2mri - convert a CTF head into a voxel coordinate
%
% MRIVertices = ctf_mri2head(HeadVertices,mri)
%
% This function converts CTF head coordinates into the CTF MRI voxel
% coordinates. The input 'HeadVertices' are head space locations; they are
% Nx3 coordinates in the CTF head space (see below). The 'mri' input
% is a struct returned by ctf_read_mri, which contains the fiducials
% and a transformation matrix between head and MRI voxel coordinates.
%
% CTF MRI volumes are 256^3 voxels, 1mm^3 each. The volume index has
% an origin at the left, anterior, superior voxel, such that:
%
% Sag increases from left to right (+X Right)
% Cor increases from anterior to posterior (+Y Posterior)
% Axi increases from superior to inferior (+Z Inferior)
%
% The head space coordinates are defined in relation to the MRI fiducial
% locations (nasion, left preauricular and right preauricular). The origin
% lies half way between the left and right preauricular points and the
% coordinate axes are given as:
%
% +X is through the nasion,
% +Y is left
% +Z is superior
%
% CTF head coordinate axes are othogonalized by first taking the cross
% product of the +X and +Y vectors (from the origin through the nasion and
% left preauricular, respectively) to find the +Z vector. Then, the Y axis
% is orthogonalized to the X/Z plane using the cross product and right hand
% rule. Hence, +Y can be slightly offset from the left preauricular.
%
% $Revision: 1.1 $ $Date: 2009-01-30 03:49:27 $
% Copyright (C) 2004 Darren L. Weber
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
% History: 04/2004, Darren.Weber_at_radiology.ucsf.edu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('in development!'); return
ver = '$Revision: 1.1 $';
fprintf('CTF_HEAD2MRI [v %s]\n',ver(11:15)); tic;
%--------------------------------------------------------------------
% Extract the coordinate transform matrices from the mri struct. These are
% designed to be left multiplied into the vertices (I don't like this
% because it means the column arrays need to be transposed below)
%T.ctfHEAD2MRI = mri.hdr.transformMatrixHead2MRI;
%T.ctfMRI2HEAD = mri.hdr.transformMatrixMRI2Head;
% OK, so I got fed up with trying to figure out how to use these
% bloody matrices and decided to calculate them myself, so
% that the entries would work in a way that I can understand!
[trans,rot] = calc_tranfer_matrix(mri);
%--------------------------------------------------------------------
fprintf('...converting from cm to mm\n');
MRIVertices = HeadVertices * 10;
%--------------------------------------------------------------------
% Pad out the HeadVertices to Nx4 matrix
Nvertices = size(MRIVertices,1);
right_column = ones( Nvertices, 1 );
MRIVertices = [ MRIVertices right_column ];
%--------------------------------------------------------------------
% Convert CTF head space coordinates into voxel indices
MRIVertices = (-1 * MRIVertices * trans) * rot;
MRIVertices = MRIVertices(:,1:3);
t = toc; fprintf('...done (%6.2f sec).\n\n',t);
return
%-------------------------------------------------------
function [trans,rot] = calc_tranfer_matrix(mri)
% these are the translations!
%ctf.mri.hdr.headOrigin_sagittal: 130
%ctf.mri.hdr.headOrigin_coronal: 123
%ctf.mri.hdr.headOrigin_axial: 152
% This is how they are calculated from the fiducials:
nas(1) = mri.hdr.HeadModel_Info.Nasion_Sag;
nas(2) = mri.hdr.HeadModel_Info.Nasion_Cor;
nas(3) = mri.hdr.HeadModel_Info.Nasion_Axi;
lpa(1) = mri.hdr.HeadModel_Info.LeftEar_Sag;
lpa(2) = mri.hdr.HeadModel_Info.LeftEar_Cor;
lpa(3) = mri.hdr.HeadModel_Info.LeftEar_Axi;
rpa(1) = mri.hdr.HeadModel_Info.RightEar_Sag;
rpa(2) = mri.hdr.HeadModel_Info.RightEar_Cor;
rpa(3) = mri.hdr.HeadModel_Info.RightEar_Axi;
headOffset_Sag = ( rpa(1) - lpa(1) ) / 2;
headOrigin_Sag = lpa(1) + headOffset_Sag;
headOffset_Cor = ( rpa(2) - lpa(2) ) / 2;
headOrigin_Cor = lpa(2) + headOffset_Cor;
headOffset_Axi = ( rpa(3) - lpa(3) ) / 2;
headOrigin_Axi = lpa(3) + headOffset_Axi;
headOrigin = [headOrigin_Sag headOrigin_Cor headOrigin_Axi];
%-------------from here to...
% calculate voxel space vectors, in slices
voxNASvector = nas - headOrigin;
voxLPAvector = lpa - headOrigin;
voxRPAvector = rpa - headOrigin;
% calculate head space vectors, in mm
headNASvector = -1 * [ voxNASvector(2), voxNASvector(1), voxNASvector(3) ];
headLPAvector = -1 * [ voxLPAvector(2), voxLPAvector(1), voxLPAvector(3) ];
headRPAvector = -1 * [ voxRPAvector(2), voxRPAvector(1), voxRPAvector(3) ];
%-------------here; is encapsulated in trans:
trans = eye(3);
trans(1,1) = 0;
trans(2,1) = 1;
trans(1,2) = 1;
trans(2,2) = 0;
trans(4,1:3) = -1 * [ headOrigin(2) headOrigin(1) headOrigin(3) ];
%headVector = -1 * ctfVox * trans;
% At this point, the voxel indices are now in the head space, in
% mm; however, the head space at this point is not orthogonal
%--------------now the head space is orthogonalized
headNASunit = headNASvector / norm(headNASvector);
headLPAunit = headLPAvector / norm(headLPAvector);
headRPAunit = headRPAvector / norm(headRPAvector);
headVERTEXunit = cross( headNASunit, headLPAunit );
% revise the LPA unit vector, so it is orthogonal to Nasion
headLPAunit = cross( headVERTEXunit, headNASunit );
% CHECK, these dot products = 0
%dot( headNASunit, headLPAunit )
%dot( headNASunit, headVERTEXunit )
%dot( headLPAunit, headVERTEXunit )
% Note that the LPA/RPA moves! This has the effect of
% rotating the coordinates in the XY plane.
rot = [ headNASunit; headLPAunit; headVERTEXunit ];
return