freesurfer2ctf_fiducial_transform

0001 
0002 clear all
0003 close all
0004 
0005 if isunix,
0006   dataPath = '/data/freesurfer/subjects/ucsf_mh/';
0007 else
0008   dataPath = 'e:\freesurfer\subjects\ucsf_mh\';
0009 end
0010 
0011 plot_fs_surf    = 0;
0012 plot_fs_voxels  = 0;
0013 plot_ctf_voxels = 0;
0014 
0015 
0016 plot_fs_fid_ras = 0;
0017 plot_ctf_fid    = 0;
0018 
0019 plot_fs2ctf_fid  = 0;
0020 plot_fs2ctf_surf = 0;
0021 
0022 
0023 % the freesurfer MRI volume (cor-???) was converted to analyze
0024 % format using the following command:
0025 %  mri_convert orig/ -oid -1 0 0 -ojd 0 1 0 -okd 0 0 1 analyze/ucsf_mh_orig_axial_las.img
0026 
0027 % This is what mri_convert calls the Analyze output matrix:
0028 % avw_mat = [-1.000   0.000   0.000   129.000;
0029 %             0.000   1.000   0.000  -129.000;
0030 %             0.000   0.000   1.000  -129.000;
0031 %             0.000   0.000   0.000     1.000;
0032 %           ];
0033 
0034 % Perhaps, by application of this matrix to the freesurfer surface,
0035 % we obtain the freesurfer vertex locations in the space of the
0036 % Analyze volume
0037 
0038 %avwFile = [dataPath,'mri/analyze/ucsf_mh_orig_axial_las.hdr'];
0039 %avw = avw_read(avwFile);
0040 
0041 
0042 %-------------------------------------------------------
0043 % Read CTF MRI volume
0044 
0045 if isunix,
0046   ctfFile = [dataPath,'mri/analyze/ucsf_mh_orig_axial_las2ctf.mri'];
0047 else
0048   ctfFile = [dataPath,'mri\analyze\ucsf_mh_orig_axial_las2ctf.mri']
0049 end
0050 
0051 % CTF MRI default volume index has an origin at the left, anterior, superior
0052 % voxel, such that Sag slices increase from left to right, Cor
0053 % slices increase from anterior to posterior, and Axi slices
0054 % increase from superior to inferior
0055 
0056 % CTF volume is 256^3 voxels, 1mm^3 each
0057 % Sag increases from left to right (+X Right)
0058 % Cor increases from anterior to posterior (+Y Posterior)
0059 % Axi increases from superior to inferior (+Z Inferior)
0060 % The origin is at the left, anterior, superior voxel
0061 % CTF volume_index is: Sag, Cor, Axi (256^3, 1mm^3 voxels)
0062 
0063 % The fiducial values can be read from the MRI volume, like this:
0064 ctf.mri = ctf_read_mri(ctfFile);
0065 headModel = ctf.mri.hdr.HeadModel_Info;
0066 
0067 
0068 %--------------------------------------------------------------------------
0069 % TESTING CODE
0070 
0071 
0072 % CTF MRI default volume index has an origin at the left, anterior, superior
0073 % voxel, such that Sag slices increase from left to right, Cor
0074 % slices increase from anterior to posterior, and Axi slices
0075 % increase from superior to inferior
0076 
0077 % CTF volume is 256^3 voxels, 1mm^3 each
0078 % Sag increases from left to right (+X Right)
0079 % Cor increases from anterior to posterior (+Y Posterior)
0080 % Axi increases from superior to inferior (+Z Inferior)
0081 % The origin is at the left, anterior, superior voxel
0082 % CTF volume_index is: Sag, Cor, Axi (256^3, 1mm^3 voxels)
0083 % CTF Sag = 256 - FreeSurfer Sag
0084 % CTF Axi = FreeSurfer Axi
0085 % CTF Cor = 256 - FreeSurfer Cor
0086 ctf.volume_index.nas = [128  35 130]; % voxels
0087 ctf.volume_index.lpa = [ 53 123 152];
0088 ctf.volume_index.rpa = [207 123 152];
0089 
0090 ctf.volume_index.mat = [ctf.volume_index.nas; ctf.volume_index.lpa; ctf.volume_index.rpa ];
0091 
0092 ctfVox = [ [ctf.volume_index.mat; 0 0 0] [ 1 1 1 0 ]' ];
0093 
0094 % eg
0095 % ctfmat =
0096 %
0097 %    128    35   130     1
0098 %     53   123   152     1
0099 %    207   123   152     1
0100 %      0     0     0     0
0101 
0102 % FreeSurfer volume_index is: Sag, Axi, Cor (256^3, 1mm^3 voxels)
0103 % FreeSurfer Sag = 256 - CTF Sag
0104 % FreeSurfer Axi = CTF Axi
0105 % FreeSurfer Cor = 256 - CTF Cor
0106 fs.volume_index.nas = [128 130 221];
0107 fs.volume_index.lpa = [203 152 133];
0108 fs.volume_index.rpa = [ 49 152 133];
0109 
0110 fs.volume_index.mat = [fs.volume_index.nas; fs.volume_index.lpa; fs.volume_index.rpa ];
0111 
0112 fsVox = [ [fs.volume_index.mat; 0 0 0] [1 1 1 0]']
0113 
0114 % eg,
0115 % fsVox =
0116 %
0117 %    128   130   221     1
0118 %    203   152   133     1
0119 %     49   152   133     1
0120 %      0     0     0     0
0121 
0122 
0123 
0124 % This is the transformation matrix from FreeSurfer voxel indices
0125 % into CTF MRI voxel coordinates.
0126 T.fs2ctf = [  -1     0     0     0  ;
0127                0     0     1     0  ;
0128                0    -1     0     0  ;
0129              256   256     0     1  ];
0130 
0131 % This is the transformation matrix from CTF voxel indices into FreeSurfer
0132 % voxel indices
0133 T.ctf2fs = [  -1     0     0     0 ;
0134                0     0    -1     0 ;
0135                0     1     0     0 ;
0136              256     0   256     1 ];
0137 
0138 % or just, T.ctf2fs = inv(T.fs2ctf);
0139 
0140 fs2ctf = fsVox  * T.fs2ctf;
0141 
0142 ctf2fs = ctfVox * T.ctf2fs;
0143 
0144 
0145 
0146 % these .sph values are not the same as ctf.mri.hdr.headOrigin_*
0147 % values; these ones correspond to
0148 % ctf.mri.hdr.HeadModel_Info.defaultSphereXYZ, which are actually
0149 % given in mm, whereas the values below in ctf.volume_xyz.sph are
0150 % in cm (taken from the MRIViewer window
0151 ctf.volume_xyz.sph = [ 0 0 5 ]; % cm, so 50 slices above left/right, in
0152                                 % the middle of left/right ear
0153 
0154 % CTF volume_xyz is: +X nasion, +Y left, +Z superior in cm
0155 % with an origin that lies half way between the left/right ears
0156 ctf.volume_xyz.nas = [ 9.07  0.00 0.00];
0157 ctf.volume_xyz.lpa = [ 0.00  7.70 0.00];
0158 ctf.volume_xyz.rpa = [ 0.00 -7.70 0.00];
0159 
0160 ctf.volume_xyz.mat = [ctf.volume_xyz.nas; ctf.volume_xyz.lpa; ctf.volume_xyz.rpa ];
0161 
0162 ctfXYZ = [ [ctf.volume_xyz.mat; 0 0 0] [ 1 1 1 0 ]' ];
0163 
0164 
0165 
0166 % END OF TESTING CODE
0167 %-------------------------------------------------------------------------
0168 
0169 
0170 
0171 
0172 
0173 
0174 
0175 
0176 
0177 
0178 % THESE ARE THE GEMS !!!
0179 
0180 %-------------------------------------------------------
0181 % Read FreeSurfer Surface
0182 
0183 if isunix,
0184   fsFile = [dataPath,'surf/rh.pial'];
0185 else
0186   fsFile = [dataPath,'surf\rh.pial'];
0187 end
0188 [fsSurf.vertices,fsSurf.faces] = freesurfer_read_surf(fsFile);
0189 
0190 % reduce the surface density for the sake of this testing
0191 % fsSurf = reducepatch(fsSurf,10000);
0192 
0193 
0194 %-------------------------------------------------------
0195 % Convert FreeSurfer Surface into FreeSurfer voxel coordinates
0196 
0197 fsVox.vertices = freesurfer_surf2voxels(fsSurf.vertices);
0198 
0199 fsVox.faces = fsSurf.faces;
0200 
0201 %-------------------------------------------------------
0202 % Convert FreeSurfer voxels into CTF voxel coordinates
0203 
0204 ctf.vertices.mri = freesurfer_surf2ctf(fsSurf.vertices);
0205 
0206 ctf.faces.mri = fsSurf.faces;
0207 
0208 ctf_write_mrishape(ctf.vertices.mri,ctf.mri);
0209 
0210 %-------------------------------------------------------
0211 % Convert CTF voxels into CTF xyz coordinates
0212 
0213 ctf.vertices.head = ctf_mri2head(ctf.vertices.mri,ctf.mri);
0214 
0215 ctf.faces.head = fsSurf.faces;
0216 
0217 ctf_write_headshape(ctf.vertices.head,ctf.mri);
0218 
0219 
0220 %-------------------------------------------------------
0221 % Write out to BrainStorm
0222 
0223 if isunix,
0224   dataPath = '/data/brainstorm_database/subjects/ucsf_mgh/';
0225 else
0226   dataPath = 'e:\brainstorm_database\subjects\ucsf_mgh\';
0227 end
0228 
0229 subjecttessFile = [dataPath,'ucsf_mgh_subjecttess'];
0230 
0231 
0232 load( [subjecttessFile,'_backup'] )
0233 figure('name','current brainstorm')
0234 patch('vertices',Vertices{1}','faces',Faces{1},'facecolor',[.75 .7 .7],'edgecolor','none');
0235 view([0,0]); light
0236 figure('name','new brainstorm')
0237 patch('vertices',ctf.vertices.head / 100,'faces',ctf.faces.head,'facecolor',[.75 .7 .7],'edgecolor','none');
0238 view([0,0]); light
0239 
0240 
0241 Vertices{1} = ctf.vertices.head' / 100;  % convert cm to meter
0242 Faces{1} = ctf.faces.head;
0243 
0244 save(subjecttessFile, 'Comment', 'Faces', 'Vertices')
0245 
0246 return
0247 
0248 
0249 
0250 
0251 
0252 
0253 
0254 
0255 
0256 
0257 
0258 %-------------------------------------------------------
0259 % Do some plots
0260 
0261 if plot_fs_surf,
0262   
0263   % To confirm that the freesurfer surface is read correctly and the
0264   % fiducials are correct, we can plot them all
0265   figure('name','FreeSurfer Surface')
0266   hold on
0267   
0268   patch('faces',fsSurf.faces,'vertices',fsSurf.vertices, ...
0269     'edgecolor','none', ...
0270     'facecolor',[.8 .7 .7],'facealpha',0.5);
0271   
0272   scatter3(fs.volume_xyz.nas(1),fs.volume_xyz.nas(2),fs.volume_xyz.nas(3),40,'r','filled');
0273   scatter3(fs.volume_xyz.lpa(1),fs.volume_xyz.lpa(2),fs.volume_xyz.lpa(3),40,'g','filled');
0274   scatter3(fs.volume_xyz.rpa(1),fs.volume_xyz.rpa(2),fs.volume_xyz.rpa(3),40,'b','filled');
0275   text(fs.volume_xyz.nas(1),fs.volume_xyz.nas(2),fs.volume_xyz.nas(3),'. nasion');
0276   text(fs.volume_xyz.lpa(1),fs.volume_xyz.lpa(2),fs.volume_xyz.lpa(3),'. left');
0277   text(fs.volume_xyz.rpa(1),fs.volume_xyz.rpa(2),fs.volume_xyz.rpa(3),'. right');
0278   
0279   set(gca,'DataAspectRatio',[1 1 1]);
0280   set(gca,'Xlim',[-128 128],'Ylim',[-128 128],'Zlim',[-128 128])
0281   rotate3d
0282   view(2)
0283   light
0284 end
0285 
0286 if plot_fs_voxels,
0287   
0288   % To confirm that the freesurfer surface is read correctly and the
0289   % fiducials are correct, we can plot them all
0290   figure('name','FreeSurfer Voxels')
0291   hold on
0292   
0293   patch('faces',fsVox.faces,'vertices',fsVox.vertices, ...
0294     'edgecolor','none', ...
0295     'facecolor',[.8 .7 .7],'facealpha',0.5);
0296   
0297   scatter3(fs.volume_index.nas(1),fs.volume_index.nas(2),fs.volume_index.nas(3),40,'r','filled');
0298   scatter3(fs.volume_index.lpa(1),fs.volume_index.lpa(2),fs.volume_index.lpa(3),40,'g','filled');
0299   scatter3(fs.volume_index.rpa(1),fs.volume_index.rpa(2),fs.volume_index.rpa(3),40,'b','filled');
0300   text(fs.volume_index.nas(1),fs.volume_index.nas(2),fs.volume_index.nas(3),'. nasion');
0301   text(fs.volume_index.lpa(1),fs.volume_index.lpa(2),fs.volume_index.lpa(3),'. left');
0302   text(fs.volume_index.rpa(1),fs.volume_index.rpa(2),fs.volume_index.rpa(3),'. right');
0303   
0304   set(gca,'DataAspectRatio',[1 1 1]);
0305   set(gca,'Xlim',[0 256],'Ylim',[0 256],'Zlim',[0 256])
0306   rotate3d
0307   view([0,0])
0308   light
0309   
0310 end
0311 
0312 if plot_ctf_voxels,
0313   
0314   % Now we can plot the CTF fiducials, which should be about 90
0315   % degree rotations of the FreeSurfer fiducials
0316   
0317   figure('name','CTF Voxels')
0318   hold on
0319   
0320   patch('faces',ctf.shape.faces,'vertices',ctf.shape.vertices, ...
0321     'edgecolor','none', ...
0322     'facecolor',[.8 .7 .7],'facealpha',0.5);
0323   
0324   scatter3(ctf.volume_index.nas(1),ctf.volume_index.nas(2),ctf.volume_index.nas(3),40,'r','filled');
0325   scatter3(ctf.volume_index.lpa(1),ctf.volume_index.lpa(2),ctf.volume_index.lpa(3),40,'g','filled');
0326   scatter3(ctf.volume_index.rpa(1),ctf.volume_index.rpa(2),ctf.volume_index.rpa(3),40,'b','filled');
0327   text(ctf.volume_index.nas(1),ctf.volume_index.nas(2),ctf.volume_index.nas(3),'. nasion');
0328   text(ctf.volume_index.lpa(1),ctf.volume_index.lpa(2),ctf.volume_index.lpa(3),'. left');
0329   text(ctf.volume_index.rpa(1),ctf.volume_index.rpa(2),ctf.volume_index.rpa(3),'. right');
0330   
0331   set(gca,'DataAspectRatio',[1 1 1]);
0332   set(gca,'Xlim',[0 256],'Ylim',[0 256],'Zlim',[0 256])
0333   rotate3d
0334   view([0,-90])
0335   light
0336 end