0001 function [stack] = fiff_read_mri(fname,read_data)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032 global FIFF;
0033 if isempty(FIFF)
0034 FIFF = fiff_define_constants();
0035 end
0036
0037 me='MNE:fiff_read_mri';
0038
0039 if nargin == 1
0040 read_data = true;
0041 elseif nargin ~= 2
0042 error(me,'Incorrect number of arguments');
0043 end
0044
0045
0046
0047 [ fid, tree ] = fiff_open(fname);
0048
0049
0050
0051
0052 mri = fiff_dir_tree_find(tree,FIFF.FIFFB_MRI);
0053 if length(mri) == 0
0054 fclose(fid);
0055 error(me,'Could not find MRI data');
0056 end
0057 mri = mri(1);
0058
0059 set = fiff_dir_tree_find(mri,FIFF.FIFFB_MRI_SET);
0060 if length(set) == 0
0061 fclose(fid);
0062 error(me,'Could not find MRI stack');
0063 end
0064 set = set(1);
0065
0066 slices = fiff_dir_tree_find(set,FIFF.FIFFB_MRI_SLICE);
0067 if length(slices) == 0
0068 fclose(fid);
0069 error(me,'Could not find MRI slices');
0070 end
0071
0072
0073
0074 tag = find_tag(mri,FIFF.FIFF_BLOCK_ID);
0075 if ~isempty(tag)
0076 stack.id = tag.data;
0077 end
0078
0079
0080
0081
0082 stack.trans = [];
0083 stack.ras_trans = [];
0084 for p = 1:set.nent
0085 kind = set.dir(p).kind;
0086 if kind == FIFF.FIFF_COORD_TRANS
0087 tag = fiff_read_tag(fid,set.dir(p).pos);
0088 if tag.data.to == FIFF.FIFFV_COORD_MRI && tag.data.from == FIFF.FIFFV_COORD_HEAD
0089 stack.trans = tag.data;
0090 elseif tag.data.to == FIFF.FIFFV_MNE_COORD_RAS && tag.data.from == FIFF.FIFFV_COORD_MRI
0091 stack.ras_trans = tag.data;
0092 end
0093 end
0094 end
0095 if isempty(stack.trans)
0096 stack.trans.from = FIFF.FIFFV_COORD_HEAD;
0097 stack.trans.to = FIFF.FIFFV_COORD_MRI;
0098 stack.trans.trans = eye(4,4);
0099 end
0100 if isempty(stack.ras_trans)
0101 stack.ras_trans.from = FIFF.FIFFV_COORD_MRI;
0102 stack.ras_trans.to = FIFF.FIFFV_MNE_COORD_RAS;
0103 stack.ras_trans.trans = eye(4,4);
0104 end
0105 stack.voxel_trans = [];
0106 stack.nslice = length(slices);
0107 if read_data
0108 fprintf(1,'\tReading slice information and pixel data.');
0109 else
0110 fprintf(1,'\tReading slice information.');
0111 end
0112 for k = 1:stack.nslice
0113 try
0114 stack.slices(k) = read_slice(slices(k));
0115 catch
0116 fclose(fid);
0117 error(me,'%s',mne_omit_first_line(lasterr));
0118 end
0119 if mod(k,50) == 0
0120 fprintf(1,'.%d.',k);
0121 end
0122 end
0123 fprintf(1,'.%d..[done]\n',k);
0124
0125 add_voxel_transform
0126
0127 fclose(fid);
0128
0129
0130 return;
0131
0132 function add_voxel_transform
0133
0134 if stack.nslice < 2
0135 fprintf(1,'\tOnly one slice in the stack. Voxel transformation will not be included\n');
0136 return;
0137 end
0138
0139
0140
0141 for k = 1:stack.nslice-1
0142 d = stack.slices(k+1).trans.trans(1:3,4)-stack.slices(k).trans.trans(1:3,4);
0143 d = sqrt(d'*d);
0144 if k == 1
0145 d0 = d;
0146 r0 = stack.slices(k).trans.trans(1:3,4);
0147 n0 = stack.slices(k).trans.trans(1:3,3);
0148 else
0149 if abs(stack.slices(k).pixel_width-stack.slices(1).pixel_width) > 1e-4 || abs(stack.slices(k).pixel_height-stack.slices(1).pixel_height) > 1e-4
0150 fprintf(1,'\tPixel sizes are not equal. Voxel transformation will not be included\n');
0151 return;
0152 end
0153 if stack.slices(k).width ~= stack.slices(1).width || stack.slices(k).height ~= stack.slices(1).height
0154 fprintf(1,'\tImage sizes are not equal. Voxel transformation will not be included\n');
0155 return;
0156 end
0157 if abs(d-d0) > 1e-4
0158 fprintf(1,'\tThe slices are not equally spaced. Voxel transformation will not be included\n');
0159 return;
0160 end
0161
0162
0163
0164 r1 = r0 + (k-1)*d0*n0;
0165 d = stack.slices(k).trans.trans(1:3,4) - r1;
0166 d = sqrt(d'*d);
0167 if abs(d) > 1e-4
0168 fprintf(1,'\tThe slices do not form a rectangular volume. Voxel transformation will not be included\n');
0169 return;
0170 end
0171 end
0172 end
0173
0174
0175
0176 voxel_trans = stack.slices(1).trans;
0177 t = eye(4);
0178 t(1,1) = stack.slices(1).pixel_width;
0179 t(2,2) = stack.slices(1).pixel_height;
0180 d = stack.slices(2).trans.trans(1:3,4)-stack.slices(1).trans.trans(1:3,4);
0181 t(3,3) = sqrt(d'*d);
0182 voxel_trans.trans = voxel_trans.trans*t;
0183 voxel_trans.from = FIFF.FIFFV_MNE_COORD_MRI_VOXEL;
0184 stack.voxel_trans = voxel_trans;
0185
0186 fprintf(1,'\tVoxel transformation added\n');
0187
0188 end
0189
0190 function [slice] = read_slice(node)
0191
0192
0193
0194 tag = find_tag(node,FIFF.FIFF_COORD_TRANS);
0195 if isempty(tag)
0196 error(me,'Could not find slice coordinate transformation');
0197 end
0198 slice.trans = tag.data;
0199 if slice.trans.from ~= FIFF.FIFFV_COORD_MRI_SLICE || ...
0200 slice.trans.to ~= FIFF.FIFFV_COORD_MRI
0201 error(me,'Illegal slice coordinate transformation');
0202 end
0203
0204
0205
0206
0207
0208
0209 slice.trans.trans(1:3,2) = -slice.trans.trans(1:3,2);
0210
0211
0212
0213 tag = find_tag(node,FIFF.FIFF_MRI_PIXEL_ENCODING);
0214 if isempty(tag)
0215 error(me,'Pixel encoding tag missing');
0216 end
0217 slice.encoding = tag.data;
0218
0219
0220
0221 tag = find_tag(node,FIFF.FIFF_MRI_PIXEL_DATA_OFFSET);
0222 if isempty(tag)
0223 slice.offset = -1;
0224 else
0225 slice.offset = tag.data;
0226 end
0227
0228
0229
0230 tag = find_tag(node,FIFF.FIFF_MRI_SOURCE_FORMAT);
0231 if isempty(tag)
0232 slice.source_format = 0;
0233 else
0234 slice.source_format = tag.data;
0235 end
0236
0237
0238
0239
0240 tag = find_tag(node,FIFF.FIFF_MRI_PIXEL_SCALE);
0241 if isempty(tag)
0242 slice.scale = 1.0;
0243 else
0244 slice.scale = tag.data;
0245 end
0246
0247
0248
0249 tag = find_tag(node,FIFF.FIFF_MRI_WIDTH);
0250 if isempty(tag)
0251 error(me,'Slice width missing');
0252 end
0253 slice.width = tag.data;
0254
0255 tag = find_tag(node,FIFF.FIFF_MRI_HEIGHT);
0256 if isempty(tag)
0257 error(me,'Slice height missing');
0258 end
0259 slice.height = tag.data;
0260
0261
0262
0263 tag = find_tag(node,FIFF.FIFF_MRI_WIDTH_M);
0264 if isempty(tag)
0265 error(me,'Pixel width missing');
0266 end
0267 slice.pixel_width = double(tag.data)/double(slice.width);
0268
0269 tag = find_tag(node,FIFF.FIFF_MRI_HEIGHT_M);
0270 if isempty(tag)
0271 error(me,'Pixel height missing');
0272 end
0273 slice.pixel_height = double(tag.data)/double(slice.height);
0274
0275
0276
0277 tag = find_tag(node,FIFF.FIFF_MRI_SOURCE_PATH);
0278 if isempty(tag)
0279 slice.offset = -1;
0280 slice.source = [];
0281
0282
0283
0284 if read_data
0285 tag = find_tag(node,FIFF.FIFF_MRI_PIXEL_DATA);
0286 if isempty(tag)
0287 error(me,'Embedded pixel data missing');
0288 end
0289 if slice.encoding == FIFF.FIFFV_MRI_PIXEL_WORD
0290 if tag.type ~= slice.encoding && tag.type ~= FIFF.FIFFT_USHORT
0291 error(me,'Embedded data is in wrong format (expected %d, got %d)',...
0292 slice.encoding,tag.type);
0293 end
0294 else
0295 if tag.type ~= slice.encoding
0296 error(me,'Embedded data is in wrong format (expected %d, got %d)',...
0297 slice.encoding,tag.type);
0298 end
0299 end
0300 if length(tag.data) ~= slice.width*slice.height
0301 error(me,'Wrong length of pixel data');
0302 end
0303
0304
0305
0306 slice.data = reshape(tag.data,slice.width,slice.height)';
0307 end
0308 else
0309 if slice.offset < 0
0310 error(me,'Offset to external file missing');
0311 end
0312 slice.source = tag.data;
0313
0314
0315
0316 if read_data
0317 pname = search_pixel_file(slice.source,fname);
0318 if isempty(pname)
0319 error(me,'Could not locate pixel file %s',slice.source);
0320 else
0321 try
0322 slice.data = read_external_pixels(pname,...
0323 slice.offset,slice.encoding,slice.width,slice.height);
0324 catch
0325 error(me,'%s',mne_omit_first_line(lasterr));
0326 end
0327 end
0328 end
0329 end
0330 end
0331
0332 function [name] = search_pixel_file(pname,sname)
0333
0334
0335
0336 if exist(pname,'file') == 2
0337 name = pname;
0338 else
0339
0340
0341
0342 a = findstr(sname,'/');
0343 if isempty(a)
0344 d = pwd;
0345 else
0346 d = sname(1:a(length(a))-1);
0347 end
0348 a = findstr(pname,'/');
0349 if ~isempty(a)
0350 pname = pname(a(length(a))+1:length(pname));
0351 end
0352 pname = sprintf('%s/../slices/%s',d,pname);
0353 if exist(pname,'file') == 2
0354 name = pname;
0355 else
0356 name = [];
0357 end
0358 end
0359 return;
0360
0361 end
0362
0363 function [pixels] = read_external_pixels(pname,offset,encoding,width,height)
0364
0365
0366
0367 if (encoding == FIFF.FIFFV_MRI_PIXEL_SWAP_WORD)
0368 sfid = fopen(pname,'rb','ieee-le');
0369 else
0370 sfid = fopen(pname,'rb','ieee-be');
0371 end
0372 if sfid < 0
0373 error(me,'Could not open pixel data file : %s',pname);
0374 end
0375 try
0376 fseek(sfid,double(offset),'bof');
0377 catch
0378 error(me,'Could not position to pixel data @ %d',offset);
0379 fclose(sfid);
0380 end
0381
0382
0383
0384 switch encoding
0385
0386 case FIFF.FIFFV_MRI_PIXEL_BYTE
0387 pixels = fread(sfid,double(width*height),'uint8=>uint8');
0388 case FIFF.FIFFV_MRI_PIXEL_WORD
0389 pixels = fread(sfid,double(width*height),'int16=>int16');
0390 case FIFF.FIFFV_MRI_PIXEL_SWAP_WORD
0391 pixels = fread(sfid,double(width*height),'int16=>int16');
0392 case FIFF.FIFFV_MRI_PIXEL_FLOAT
0393 pixels = fread(sfid,double(width*height),'single=>double');
0394 otherwise
0395 fclose(sfid);
0396 error(me,'Unknown pixel encoding : %d',encoding);
0397 end
0398 fclose(sfid);
0399
0400
0401
0402 pixels = reshape(pixels,width,height)';
0403 end
0404
0405 function [tag] = find_tag(node,findkind)
0406
0407 for p = 1:node.nent
0408 kind = node.dir(p).kind;
0409 pos = node.dir(p).pos;
0410 if kind == findkind
0411 tag = fiff_read_tag(fid,pos);
0412 return;
0413 end
0414 end
0415 tag = [];
0416 return
0417
0418 end
0419
0420 end
0421