0001 function d = SurfStatReadVol1( file, Z, T );
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 if ~isstr(file)
0014 d = file;
0015 if isfield(d,'data') & nargin==3
0016 if length(size(d.data))==4 & Z~=0 & T~=0
0017 d.data=file.data(:,:,Z,T);
0018 end
0019 if length(size(d.data))==3 & length(Z)>1 & length(T)==1
0020 d.data=file.data(:,:,Z);
0021 end
0022 if length(size(d.data))==3 & length(Z)==1 & length(T)>1
0023 d.data=file.data(:,:,T);
0024 end
0025 if length(size(d.data))==3 & length(Z)==1 & length(T)==1 & Z>0 & T>0
0026 d.data=file.data(:,:,Z);
0027 end
0028 end
0029 return
0030 end
0031
0032 file = deblank(file);
0033
0034 [path,name,ext]=fileparts(file);
0035 if strcmp(ext,'.gz')
0036 if isunix
0037 unix(['gunzip ' file]);
0038 else
0039 ['Can''t gunzip on non-unix system']
0040 return
0041 end
0042 [path,name,ext]=fileparts(name)
0043 end
0044
0045 switch lower(ext)
0046 case '.mnc'
0047 if nargin == 3
0048 d = fmris_read_minc(file,Z,T);
0049 else
0050 d = fmris_read_minc(file);
0051 end
0052 case '.img'
0053 if nargin == 3
0054 d = fmris_read_analyze(file,Z,T);
0055 else
0056 d = fmris_read_analyze(file);
0057 end
0058 if isfield(d,'data') d.data(isnan(d.data))=0; end;
0059 d.origin=-d.origin.*d.vox(1:3);
0060 case '.brik'
0061 if nargin == 3
0062 d = fmris_read_afni(file,Z,T);
0063 else
0064 d = fmris_read_afni(file);
0065 end
0066 case '.nii'
0067 if nargin == 3
0068 d = fmris_read_nifti(file,Z,T);
0069 else
0070 d = fmris_read_nifti(file);
0071 end
0072 otherwise
0073 ['Unknown file extension']
0074 end
0075
0076 d.parent_file=file;
0077
0078 return
0079 end
0080
0081
0082
0083
0084
0085 function [d]=fmris_read_afni(file,Z,T)
0086
0087 [err, Info]=BrikInfo(file);
0088
0089 d.dim = [Info.DATASET_DIMENSIONS(1:3) Info.DATASET_RANK(2)];
0090
0091 if nargin == 1
0092 Z = 1:d.dim(3);
0093 T = 1:d.dim(4);
0094 end
0095
0096 if (T(1)~=0)&(Z(1)~=0)
0097 Opt.Slices=Z;
0098 Opt.Frames=T;
0099 [err, d.data, Info, ErrMessage] = BrikLoad(file, Opt);
0100 d.calib = [min(min(min(min(d.data)))) max(max(max(max(d.data))))];
0101 end
0102
0103 d.vox = Info.DELTA;
0104 d.vox_units = '';
0105 d.vox_offset = 0;
0106 d.precision = '';
0107 d.calib_units = '';
0108 d.origin = Info.ORIGIN;
0109 d.descrip = '';
0110
0111 if isfield(Info,'WORSLEY_DF')
0112 df=Info.WORSLEY_DF;
0113 else
0114 df=[];
0115 end
0116 if ~isempty(df)
0117 d.df=df;
0118 end
0119
0120 if isfield(Info,'WORSLEY_NCONJ')
0121 nconj=Info.WORSLEY_NCONJ;
0122 else
0123 nconj=[];
0124 end
0125 if ~isempty(nconj)
0126 d.nconj=nconj;
0127 end
0128
0129 if isfield(Info,'WORSLEY_FWHM')
0130 fwhm=Info.WORSLEY_FWHM;
0131 else
0132 fwhm=[];
0133 end
0134 if ~isempty(fwhm)
0135 d.fwhm=fwhm;
0136 end
0137
0138 return;
0139 end
0140
0141
0142
0143
0144
0145 function d=fmris_read_analyze(file,Z,T);
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162 machineformat=getmachineformat(file);
0163
0164 if ~isempty(machineformat)
0165
0166
0167 d=Read_Analyze_Hdr(file,machineformat);
0168 d.file_name=file;
0169
0170
0171 if strcmp(d.precision,'Unknown')
0172 d.precision=getprecision(deblank(file),machineformat,d.dim,d.global);
0173 end
0174
0175
0176
0177 if ~isempty(d) & ~strcmp(d.precision,'Unknown')
0178
0179 if nargin==1 | strcmp(d.precision,'uint1')
0180
0181
0182 fid = fopen(d.file_name,'r',machineformat);
0183 if fid > -1
0184 d.data=d.scale*reshape(fread(fid,prod(d.dim),d.precision),d.dim(1),d.dim(2),d.dim(3),d.dim(4));
0185 fclose(fid);
0186 else
0187 errordlg('Check Image File: Existence, Permissions ?','Read Error');
0188 end;
0189
0190 if nargin==3
0191 if all(Z>0)&all(Z<=d.dim(3))&all(T>0)&all(T<=d.dim(4))
0192 d.data=d.data(:,:,Z,T);
0193 d.Z=Z;
0194 d.T=T;
0195 else
0196 errordlg('Incompatible Matrix Identifiers !','Read Error');
0197 end
0198 end
0199
0200 elseif nargin==3
0201
0202 if (T(1)~=0)|(Z(1)~=0)
0203
0204 if all(Z>0)&all(Z<=d.dim(3))&all(T>0)&all(T<=d.dim(4))
0205
0206 fid = fopen(d.file_name,'r',machineformat);
0207 if fid > -1
0208 d.data=zeros(d.dim(1),d.dim(2),length(Z),length(T));
0209 for t=1:length(T)
0210 for z=1:length(Z)
0211 status=fseek(fid,d.hdr.byte*((T(t)-1)*prod(d.dim(1:3))+(Z(z)-1)*prod(d.dim(1:2))),'bof');
0212 d.data(:,:,z,t)=d.scale*fread(fid,[d.dim(1) d.dim(2)],d.precision);
0213 end
0214 end
0215 d.Z=Z;
0216 d.T=T;
0217 fclose(fid);
0218 else
0219 errordlg('Check Image File: Existence, Permissions ?','Read Error');
0220 end;
0221
0222 else
0223 errordlg('Incompatible Matrix Identifiers !','Read Error');
0224 end;
0225
0226 end;
0227 else
0228 errordlg('Unusual Number of Arguments','Read Error');
0229 end;
0230 else
0231 if strcmp(d.precision,'Unknown');
0232 errordlg('Unknown Data Type (Precision?)','Read Error');
0233 end
0234 end;
0235 else
0236 errordlg('Unknown Machine Format','Read Error');
0237 end
0238
0239
0240
0241 if d.vox(3)==0
0242 d.vox(3)=6;
0243 end
0244
0245 return;
0246 end
0247
0248
0249
0250 function machineformat=getmachineformat(file);
0251
0252
0253
0254
0255 machineformat=[];
0256 for mf='nlbdgcas'
0257 fid = fopen([file(1:(length(file)-3)) 'hdr'],'r',mf);
0258 if fid > -1
0259 fseek(fid,40,'bof');
0260 if any(fread(fid,1,'int16')==1:4)
0261 machineformat=mf;
0262 fclose(fid);
0263 break
0264 else
0265 fclose(fid);
0266 end
0267 else
0268 errordlg('Check Header File: Existence, Permissions ?','Read Error');
0269 break
0270 end
0271 end
0272
0273 return
0274 end
0275
0276
0277
0278 function precision=getprecision(file,machineformat,dim,range);
0279
0280
0281
0282
0283 precisions=['int8 '
0284 'int16 '
0285 'int32 '
0286 'int64 '
0287 'uint8 '
0288 'uint16 '
0289 'uint32 '
0290 'uint64 '
0291 'single '
0292 'float32'
0293 'double '
0294 'float64' ];
0295 nbytes=[1 2 4 8 1 2 4 8 4 4 8 8];
0296 middle_vol=dim(1)*dim(2)*floor(dim(3)/2)+dim(1)*round(dim(2)/2)+round(dim(1)/2);
0297 h=dir(file);
0298 n=ceil(h.bytes/prod(dim));
0299
0300 fid = fopen(file,'r',machineformat);
0301 if fid > -1
0302 for i=1:size(precisions,1)
0303 if nbytes(i)==n
0304 status=fseek(fid,middle_vol*n,'bof');
0305 if status==0
0306 precision=deblank(precisions(i,:));
0307 x=fread(fid,10,precision);
0308 if all(range(1)<=x) & all(x<=range(2))
0309 return
0310 end
0311 end
0312 end
0313 end
0314 end
0315 errordlg('Check Header File: Existence, Permissions ?','Read Error');
0316 precision='Unknown';
0317
0318 return
0319 end
0320
0321
0322
0323 function d=Read_Analyze_Hdr(file,machineformat);
0324
0325
0326
0327
0328 fid = fopen([file(1:(length(file)-3)) 'hdr'],'r',machineformat);
0329
0330 if fid > -1
0331
0332
0333
0334 fseek(fid,0,'bof');
0335
0336 d.hdr.sizeof_hdr = fread(fid,1,'int32');
0337 d.hdr.data_type = deblank(setstr(fread(fid,10,'char'))');
0338 d.hdr.db_name = deblank(setstr(fread(fid,18,'char'))');
0339 d.hdr.extents = fread(fid,1,'int32');
0340 d.hdr.session_error = fread(fid,1,'int16');
0341 d.hdr.regular = deblank(setstr(fread(fid,1,'char'))');
0342 d.hdr.hkey_un0 = deblank(setstr(fread(fid,1,'char'))');
0343
0344
0345
0346
0347
0348 fseek(fid,40,'bof');
0349
0350 d.hdr.dim = fread(fid,8,'int16');
0351
0352 d.hdr.vox_units = deblank(setstr(fread(fid,4,'char'))');
0353 d.hdr.cal_units = deblank(setstr(fread(fid,8,'char'))');
0354 d.hdr.unused1 = fread(fid,1,'int16');
0355 d.hdr.datatype = fread(fid,1,'int16');
0356 d.hdr.bitpix = fread(fid,1,'int16');
0357 d.hdr.dim_un0 = fread(fid,1,'int16');
0358 d.hdr.pixdim = fread(fid,8,'float');
0359 d.hdr.vox_offset = fread(fid,1,'float');
0360 d.hdr.funused1 = fread(fid,1,'float');
0361 d.hdr.funused2 = fread(fid,1,'float');
0362 d.hdr.funused3 = fread(fid,1,'float');
0363 d.hdr.cal_max = fread(fid,1,'float');
0364 d.hdr.cal_min = fread(fid,1,'float');
0365 d.hdr.compressed = fread(fid,1,'int32');
0366 d.hdr.verified = fread(fid,1,'int32');
0367 d.hdr.glmax = fread(fid,1,'int32');
0368 d.hdr.glmin = fread(fid,1,'int32');
0369
0370
0371
0372 fseek(fid,148,'bof');
0373
0374 d.hdr.descrip = deblank(setstr(fread(fid,80,'char'))');
0375 d.hdr.aux_file = deblank(setstr(fread(fid,24,'char'))');
0376 d.hdr.orient = fread(fid,1,'char');
0377 d.hdr.origin = fread(fid,5,'uint16');
0378 d.hdr.generated = deblank(setstr(fread(fid,10,'char'))');
0379 d.hdr.scannum = deblank(setstr(fread(fid,10,'char'))');
0380 d.hdr.patient_id = deblank(setstr(fread(fid,10,'char'))');
0381 d.hdr.exp_date = deblank(setstr(fread(fid,10,'char'))');
0382 d.hdr.exp_time = deblank(setstr(fread(fid,10,'char'))');
0383 d.hdr.hist_un0 = deblank(setstr(fread(fid,3,'char'))');
0384 d.hdr.views = fread(fid,1,'int32');
0385 d.hdr.vols_added = fread(fid,1,'int32');
0386 d.hdr.start_field = fread(fid,1,'int32');
0387 d.hdr.field_skip = fread(fid,1,'int32');
0388 d.hdr.omax = fread(fid,1,'int32');
0389 d.hdr.omin = fread(fid,1,'int32');
0390 d.hdr.smax = fread(fid,1,'int32');
0391 d.hdr.smin = fread(fid,1,'int32');
0392
0393 fclose(fid);
0394
0395
0396
0397
0398 d.dim = d.hdr.dim(2:5)';
0399 vox = d.hdr.pixdim(2:5)';
0400 if vox(4)==0
0401 vox(4)=[];
0402 end
0403 d.vox = vox;
0404 d.vox_units = d.hdr.vox_units;
0405 d.vox_offset = d.hdr.vox_offset;
0406 scale = d.hdr.funused1;
0407 d.scale = ~scale + scale;
0408 d.global = [d.hdr.glmin d.hdr.glmax];
0409 d.calib = [d.hdr.cal_min d.hdr.cal_max];
0410 d.calib_units = d.hdr.cal_units;
0411 d.origin = d.hdr.origin(1:3)';
0412 d.descrip = d.hdr.descrip(1:max(find(d.hdr.descrip)));
0413
0414 if ~isnan(str2double(d.hdr.scannum)); d.df(1) =str2double(d.hdr.scannum); end
0415 if ~isnan(str2double(d.hdr.patient_id)); d.df(2) =str2double(d.hdr.patient_id); end
0416 if ~isnan(str2double(d.hdr.exp_date)); d.fwhm(1)=str2double(d.hdr.exp_date); end
0417 if ~isnan(str2double(d.hdr.exp_time)); d.fwhm(2)=str2double(d.hdr.exp_time); end
0418 if ~isnan(str2double(d.hdr.hist_un0)); d.nconj =str2double(d.hdr.hist_un0); end
0419
0420 switch d.hdr.datatype
0421 case 1
0422 d.precision = 'uint1';
0423 d.hdr.byte = 0;
0424 case 2
0425 d.precision = 'uint8';
0426 d.hdr.byte = 1;
0427 case 4
0428 d.precision = 'int16';
0429 d.hdr.byte = 2;
0430 case 8
0431 d.precision = 'int32';
0432 d.hdr.byte = 4;
0433 case 16
0434 d.precision = 'float';
0435 d.hdr.byte = 4;
0436 case 64
0437 d.precision = 'double';
0438 d.hdr.byte = 8;
0439 otherwise
0440 d.precision = 'Unknown';
0441 d.hdr.byte = 0;
0442 end
0443
0444 else
0445 d=[];
0446 errordlg('Check Header File: Existence, Permissions ?','Read Error');
0447 end
0448
0449 return
0450 end
0451
0452
0453
0454
0455 function [d]=fmris_read_minc(file,Z,T)
0456
0457 fid=fopen(file);
0458 file=fopen(fid);
0459 fclose(fid);
0460
0461 d.file_name = file;
0462
0463 d.dim = fliplr(miinquire(file,'imagesize')');
0464 d.dim = d.dim + (d.dim == 0);
0465
0466 if nargin == 1
0467 Z = 1:d.dim(3);
0468 T = 1:d.dim(4);
0469 end
0470
0471 if (T~=0)&(Z~=0)
0472 if d.dim(4)==1
0473 images = mireadimages(file,Z-1);
0474 else
0475 n=length(T);
0476 images=[];
0477 for i=0:floor((n-1)/256)
0478 frames=T((i*256+1):min((i+1)*256,n));
0479 images = [images mireadimages(file,Z-1,frames-1)];
0480 end
0481 end
0482 d.data = reshape(images,[d.dim(1:2) length(Z) length(T)]);
0483 d.calib = [min(min(min(min(d.data)))) max(max(max(max(d.data))))];
0484 end
0485
0486 [x_vox, y_vox, z_vox] = miinquire(file,'attvalue','xspace','step','attvalue','yspace','step','attvalue','zspace','step');
0487 d.vox = [x_vox y_vox z_vox];
0488 d.vox_units = '';
0489 d.vox_offset = 0;
0490 d.precision = '';
0491 d.calib_units = '';
0492 [x_origin, y_origin, z_origin] = miinquire(file,'attvalue','xspace','start','attvalue','yspace','start','attvalue','zspace','start');
0493 d.origin = [x_origin y_origin z_origin];
0494 d.descrip = '';
0495
0496 df=miinquire(file, 'attvalue', 'df', 'df');
0497 if ~isempty(df)
0498 d.df=df;
0499 end
0500
0501 nconj=miinquire(file, 'attvalue', 'nconj', 'nconj');
0502 if ~isempty(nconj)
0503 d.nconj=nconj;
0504 end
0505
0506 fwhm=miinquire(file, 'attvalue', 'fwhm', 'fwhm');
0507 if ~isempty(fwhm)
0508 d.fwhm=fwhm;
0509 end
0510
0511 return;
0512 end
0513
0514
0515
0516
0517 function [d]=fmris_read_nifti(file,Z,T)
0518
0519 d.file_name=file;
0520 machineformats='nlbdgcas';
0521 for i=1:length(machineformats)
0522 fid=fopen(file,'r',machineformats(i));
0523 if fid == -1;
0524 fprintf('Invalid file %s',file)
0525 end
0526 sizeof_hdr=fread(fid,1,'int');
0527 if sizeof_hdr==348;
0528
0529 break;
0530 else
0531
0532
0533 fclose(fid);
0534 end
0535 end
0536 data_type=fread(fid,10,'char');
0537 db_name=fread(fid,18,'char');
0538 extents=fread(fid,1,'int');
0539 session_error=fread(fid,1,'short');
0540 regular=char(fread(fid,1,'char')');
0541 dim_info=char(fread(fid,1,'char')');
0542 dim=fread(fid,8,'short');
0543 intent_p =fread(fid,2,'float');
0544 intent_q =fread(fid,2,'uint16');
0545 intent_code =fread(fid,1,'short');
0546 datatype=fread(fid,1,'short');
0547 bitpix=fread(fid,1,'short');
0548 slice_start=fread(fid,1,'short');
0549 pixdim=fread(fid,8,'float');
0550 vox_offset=fread(fid,1,'float');
0551 scl_slope =fread(fid,1,'float');
0552 scl_inter =fread(fid,1,'float');
0553 slice_end=fread(fid,1,'short');
0554 slice_code =char(fread(fid,1,'char')');
0555 xyzt_units =char(fread(fid,1,'char')');
0556 cal_max=fread(fid,1,'float');
0557 cal_min=fread(fid,1,'float');
0558 slice_duration=fread(fid,1,'float');
0559 toffset=fread(fid,1,'float');
0560 glmax=fread(fid,1,'int');
0561 glmin=fread(fid,1,'int');
0562 descrip=char(fread(fid,80,'char')');
0563 aux_file=char(fread(fid,24,'char')');
0564 qform_code =fread(fid,1,'short');
0565 sform_code =fread(fid,1,'short');
0566 quatern_b =fread(fid,1,'float');
0567 quatern_c =fread(fid,1,'float');
0568 quatern_d =fread(fid,1,'float');
0569 qoffset_x =fread(fid,1,'float');
0570 qoffset_y =fread(fid,1,'float');
0571 qoffset_z =fread(fid,1,'float');
0572 srow_x =fread(fid,4,'float');
0573 srow_y =fread(fid,4,'float');
0574 srow_z =fread(fid,4,'float');
0575 intent_name=char(fread(fid,16,'char')');
0576 magic =char(fread(fid,4,'char')');
0577
0578 d.machineformat=machineformats(i);
0579 d.dim=ones(1,4);
0580 if dim(1)==5
0581 if dim(5)==1
0582 dim=[4; dim(2:4); dim(6); zeros(3,1)]
0583 else
0584 dim
0585 fclose(fid);
0586 return
0587 end
0588 end
0589 d.dim(1:dim(1))=dim((1:dim(1))+1);
0590 d.vox=zeros(1,3);
0591 d.vox(1:dim(1))=pixdim((1:dim(1))+1);
0592
0593 d.vox_offset = vox_offset;
0594 d.scale = scl_slope;
0595 d.intercept = scl_inter;
0596 d.global = [glmin glmax];
0597 d.calib = [cal_min cal_max];
0598 if qform_code>0;
0599 d.origin = [qoffset_x qoffset_y qoffset_z];
0600 else
0601 d.origin = [srow_x(4) srow_y(4) srow_z(4)];
0602 end
0603 d.descrip = descrip;
0604 d.parent_file=file;
0605
0606 if intent_p(1)>0; d.df(1)=intent_p(1); end;
0607 if intent_p(2)>0; d.df(2)=intent_p(2); end;
0608
0609 intent_q=intent_q/2^16*100;
0610 if intent_q(1)>0; d.fwhm(1)=intent_q(1); end;
0611 if intent_q(2)>0; d.fwhm(2)=intent_q(2); end;
0612
0613 if nargin<2
0614 Z=1:d.dim(3);
0615 T=1:d.dim(4);
0616 end
0617
0618 if Z(1)==0 & T(1)==0
0619 fclose(fid);
0620 return
0621 end
0622
0623 types=lower(['UINT8 ';'INT16 ';'INT32 ';'FLOAT32';'FLOAT64']);
0624 codes=[2; 4; 8; 16; 64];
0625 d.precision=[];
0626 for i=1:5
0627 if datatype==codes(i)
0628 d.precision=deblank(types(i,:));
0629 end
0630 end
0631 if isempty(d.precision)
0632 unknown_datatype=datatype
0633 fclose(fid);
0634 return
0635 end
0636
0637
0638 d.byte=4;
0639 if all(Z>0)&all(Z<=d.dim(3))&all(T>0)&all(T<=d.dim(4))
0640 d.data=zeros(d.dim(1),d.dim(2),length(Z),length(T));
0641 for t=1:length(T)
0642 for z=1:length(Z)
0643 position=d.byte*((T(t)-1)*prod(d.dim(1:3))+(Z(z)-1)*prod(d.dim(1:2)))+vox_offset;
0644 status=fseek(fid,position,'bof');
0645 d.data(:,:,z,t)=fread(fid,[d.dim(1) d.dim(2)],d.precision);
0646 end
0647 end
0648 fclose(fid);
0649 else
0650 Z
0651 T
0652 fclose(fid);
0653 return
0654 end
0655
0656 if scl_slope~=0
0657 d.data=d.data*scl_slope+scl_inter;
0658 end
0659
0660
0661 return;
0662 end
0663