00001 function [err, V, Info, ErrMessage] = BrikLoad (BrikName, param1, param2)
00002 %
00003 % OLD USAGE for backward compatibility, see New Usage
00004 %
00005 % [err, V, Info, ErrMessage] = BrikLoad (BrikName, [form], [MachineFormat])
00006 %
00007 %Purpose:
00008 % loads an AFNI brik into V
00009 %
00010 %
00011 %Input Parameters:
00012 % BrikName, name of the brik
00013 % form, 'vector' or 'matrix' this is optional, the default is matrix
00014 % see .Format help under NEW USAGE
00015 % MachineFormat is a string such as 'native' or 'ieee-le' (LSB_FIRST) or 'ieee-be' (MSB_FIRST)
00016 % see .MachineFormat help under NEW USAGE
00017 %
00018 %
00019 % NEW USAGE
00020 %
00021 % [err, V, Info, ErrMessage] = BrikLoad (BrikName, [Opt])
00022 %
00023 %Purpose:
00024 % loads an AFNI brik or a 1D file into V
00025 %
00026 %
00027 %Input Parameters:
00028 % BrikName, name of the brik
00029 % Opt is an options format with the following optional fields
00030 % .Format, 'vector' or 'matrix', the default is matrix
00031 % This determines how the brick data is returned in V.
00032 % If you choose 'vector' then a N x M x K volume is stored in a
00033 % (N * M * K) column vector. If the brick has multiple volumes (4-D)
00034 % then an N x M x K x J brick is stored in an (N * M * K) x J matrix.
00035 % If you use 'vector' option you can change V to matrix format by using
00036 % M = reshape(V, Info.DATASET_DIMENSIONS(1), Info.DATASET_DIMENSIONS(2),...
00037 % Info.DATASET_DIMENSIONS(3), Info.DATASET_RANK(2));
00038 %
00039 % Note that indexing in matlab is different than in AFNI. Matlab starts indexing at
00040 % 1 while AFNI starts with 0. So voxel (i,j,k) in AFNI corresponds to voxel (i+1,j+1,k+1)
00041 % in matlab. (see example below).
00042 %
00043 % .MachineFormat is a string such as 'native' or 'ieee-le' (LSB_FIRST) or 'ieee-be' (MSB_FIRST)
00044 % default is whatever is specified in the .HEAD file.
00045 % If nothing is specified in the .HEAD file, MSB_FIRST or 'ieee-be' is assumed
00046 % since most files created with the old versions of AFNI were on the SGIs .
00047 % You must specify the parameter Opt.Format, to specify Opt.MachineFormat and override the defaults
00048 % see help fopen for more info
00049 %
00050 % .Scale 0/1 if 1, then the scaling factor is applied to the values read from the brik
00051 % default is 1
00052 %
00053 % .Slices: vector of slices, 1 based. Default is all slices.
00054 % Read the set of slices specified in .Slices (added for FMRISTAT)
00055 % .SliceSize_1D: If you are reading 1D files in chunks, specify the
00056 % number of rows you want read in at any one time.
00057 % If chunk is 1000 and .Slices = 1 then you would
00058 % get values from rows 1 (the first row) to row 1000.
00059 % When .Slices = 2, you would get rows 1001 ... 2000 .
00060 % If the number of rows in your dataset is not an
00061 % integer multiple of SliceSize_1D the last read
00062 % will return the left over rows.
00063 % .Frames: vector of frames, 1 based. Default is all frames.
00064 % Read the sub-bricks specified in .Frames (added for FMRISTAT)
00065 % .method: method option for Read_1D if you are using 1D files.
00066 % see Read_1D -help for more info
00067 % WARNING: If you use .Scale = 0 , you may get bad time series if you have multiple subbriks !
00068 % Each subbrik can have a different scaling factor
00069 %
00070 %Output Parameters:
00071 % err : 0 No Problem
00072 % : 1 Mucho Problems
00073 % V : the brik data in vector or matrix (X,Y,Z,t) form
00074 % Info : is a structure containing some Brik infor (from .HEAD file)
00075 % ErrMessage: An error or warning message
00076 %
00077 %Key Terms:
00078 %
00079 %More Info :
00080 % How to read a brick and display a time series:
00081 % %read the 3D+time brick
00082 % [err, V, Info, ErrMessage] = BrikLoad ('ARzs_CW_r5+orig');
00083 % %plot the time course of voxel (29, 33, 3)
00084 % plot (squeeze(V(30 , 34, 4, :))); %squeeze is used to turn the 1x1x1x160 matrix into a 160x1 vector for the plot function
00085 %
00086 %
00087 %
00088 % see also BrikInfo, Info_1D, WriteBrik
00089 %
00090 % The core for this function (fread and reshape) is based on
00091 % Timothy M. Ellmore's (Laboratory of Brain and Cognition, NIMH)
00092 % function ReadBRIK
00093 %
00094 % .Slices and .Frames options were added by Dr. Keith Worsley for FMRISTAT
00095 % http:
00096 % http:
00097 %
00098 % Author : Ziad Saad Mon Oct 18 14:47:32 CDT 1999
00099 % Biomedical Engineering, Marquette University
00100 % Biophysics Research institute, Medical College of Wisconsin
00101 %
00102 % Contact: ziad@nih.gov
00103
00104
00105 %Define the function name for easy referencing
00106 FuncName = 'BrikLoad';
00107
00108 %Debug Flag
00109 DBG = 1;
00110
00111 %initailize return variables
00112 err = 1;
00113 V = [];
00114 Info = [];
00115 ErrMessage = '';
00116
00117 switch nargin,
00118 case 1,
00119 DoVect = 0;
00120 Opt.MachineFormat = '';
00121 Opt.Format = 'matrix';
00122 Opt.Scale = 1;
00123 Opt.Slices = [];
00124 Opt.Frames = [];
00125 Opt.FileFormat = '';
00126 case 2,
00127 if (~isstruct(param1)),
00128 Opt.Format = param1;
00129 Opt.MachineFormat = '';
00130 Opt.Scale = 1;
00131 Opt.Slices = [];
00132 Opt.Frames = [];
00133 Opt.FileFormat = '';
00134 else
00135 Opt = param1;
00136 if (~isfield(Opt,'FileFormat')), Opt.FileFormat = ''; end
00137 end
00138 case 3,
00139 Opt.Format = param1;
00140 Opt.MachineFormat = param2;
00141 Opt.Scale = 1;
00142 Opt.Slices = [];
00143 Opt.Frames = [];
00144 Opt.FileFormat = '';
00145 otherwise,
00146 err= ErrEval(FuncName,'Err_Bad option');
00147 return;
00148 end
00149
00150 % Are we dealing with a .1D file ?
00151 is1D = 0;
00152 if (isempty(Opt.FileFormat)), % try to guess
00153 [St, xtr] = Remove1DExtension(BrikName);
00154 if (~isempty(xtr)),
00155 is1D = 1;
00156 end
00157 elseif (strcmp(Opt.FileFormat, '1D') | strcmp(Opt.FileFormat, '1d')),
00158 is1D = 1;
00159 end
00160
00161 if (is1D), % 1D land
00162 V = []; Info = []; ErrMessage = '';
00163 Opt.verb = 1;
00164 if (~isfield(Opt, 'method') | isempty(Opt.method)), Opt.method = 0; end
00165 if (isfield(Opt,'Slices') & ~isempty(Opt.Slices)),
00166 if (length(Opt.Slices) ~= 1),
00167 ErrMessage = sprintf ('%s: Opt.Slices can only be used to specify one slice at a time for 1D format', FuncName, BrikName);
00168 err = ErrEval(FuncName,'Err_1D Bad .Slices option');
00169 return;
00170 end
00171 if (~isfield(Opt, 'SliceSize_1D') | isempty(Opt.SliceSize_1D)),
00172 ErrMessage = sprintf ('%s: SliceSize_1D must be specified with Slices option for 1D files', FuncName);
00173 err = ErrEval(FuncName,'Err_1D Bad .SliceSize_1D option');
00174 return;
00175 end
00176 Opt.chunk_size = Opt.SliceSize_1D;
00177 Opt.chunk_index = Opt.Slices - 1;
00178 if (isfield(Opt,'Frames') & ~isempty(Opt.Frames)), Opt.col_index = Opt.Frames - 1; end
00179 end
00180
00181 [err, V, Info] = Read_1D(BrikName, Opt);
00182 if (err),
00183 ErrMessage = sprintf ('%s: Failed to read %s file', FuncName, BrikName);
00184 err = ErrEval(FuncName,'Err_1D file could not be read');
00185 return;
00186 end
00187 return;
00188 end
00189
00190 %assume you have a brik format and proceed as usual
00191 %make sure Opt fields are set OK. That's done for the new usage format
00192 % if (~isfield(Opt,'') | isempty(Opt.)), Opt. = ; end
00193 if (~isfield(Opt,'Format') | isempty(Opt.Format)), Opt.Format = 'matrix'; end
00194 if (~isfield(Opt,'MachineFormat') | isempty(Opt.MachineFormat)), Opt.MachineFormat = ''; end
00195 if (~isfield(Opt,'Scale') | isempty(Opt.Scale)), Opt.Scale = 1; end %you can change the default for the new usage here
00196 if (~isfield(Opt,'Slices') | isempty(Opt.Slices)), Opt.Slices = []; end
00197 if (~isfield(Opt,'Frames') | isempty(Opt.Frames)), Opt.Frames = []; end
00198 if (~isfield(Opt,'FileFormat') | isempty(Opt.FileFormat)), Opt.FileFormat = ''; end
00199
00200
00201 %set the format
00202 if (eq_str(Opt.Format,'vector')),
00203 DoVect = 1;
00204 elseif(eq_str(Opt.Format,'matrix')),
00205 DoVect = 0;
00206 end
00207
00208 %Fix the name of the brik
00209 %make sure there are no .BRIK or .HEAD
00210 vtmp = findstr(BrikName,'.BRIK');
00211 if (~isempty(vtmp)), %remove .BRIK
00212 BrikName = BrikName(1:vtmp(1)-1);
00213 end
00214
00215 vtmp = findstr(BrikName,'.HEAD');
00216 if (~isempty(vtmp)), %remove .HEAD
00217 BrikName = BrikName(1:vtmp(1)-1);
00218 end
00219
00220
00221 %Now make sure .HEAD and .BRIK are present
00222 sHead = sprintf('%s.HEAD', BrikName);
00223 sBRIK = sprintf('%s.BRIK', BrikName);
00224 if (~filexist(sHead)),
00225 ErrMessage = sprintf ('%s: %s not found', FuncName, sHead);
00226 err = ErrEval(FuncName,'Err_HEAD file not found');
00227 return;
00228 end
00229 if (~filexist(sBRIK)),
00230 ErrMessage = sprintf ('%s: %s not found', FuncName, sBRIK);
00231 err = ErrEval(FuncName,'Err_BRIK file not found');
00232 return;
00233 end
00234
00235
00236 %get Brik info and setup for slice and frame selectors (for FMRISTAT)
00237 [err, Info] = BrikInfo(BrikName);
00238
00239 allslices=1:Info.DATASET_DIMENSIONS(3);
00240 if (~isfield(Opt,'Slices') | isempty(Opt.Slices)),
00241 Opt.Slices = allslices;
00242 end
00243 isallslices=all(ismember(allslices,Opt.Slices));
00244 allframes=1:Info.DATASET_RANK(2);
00245 if (~isfield(Opt,'Frames') | isempty(Opt.Frames)),
00246 Opt.Frames = allframes;
00247 end
00248 isallframes=all(ismember(allframes,Opt.Frames));
00249
00250 if (~strcmp(Info.ByteOrder,'unspecified')),
00251 %found byte order specs, make sure it does not conflict with the user option
00252 if (~isempty(Opt.MachineFormat)),
00253 %user specified a format
00254 if (~strcmp(Info.ByteOrder,Opt.MachineFormat)),
00255 %clash, warn
00256 ErrEval(FuncName,'Wrn_Machine format specified conflicts with .HEAD info. Proceeding with fread ...');
00257 end
00258 else
00259 Opt.MachineFormat = Info.ByteOrder;
00260 end
00261 else
00262 %did not find ByteOrder in .HEAD, use user option or pick default
00263 if (isempty(Opt.MachineFormat)),
00264 Opt.MachineFormat = 'ieee-be';
00265 ErrEval(FuncName,'Wrn_No Machine Format was specified by user or found in .HEAD file. Using ieee-be (MSB_FIRST)');
00266 end
00267 end
00268
00269 %figure out storage type
00270 itype = unique(Info.BRICK_TYPES);
00271 if (length(itype) > 1),
00272 err = 1; ErrMessage = sprintf('Error %s: Not set up to read sub-bricks of multiple sub-types', FuncName); errordlg(ErrMessage);
00273 return;
00274 end
00275 switch itype,
00276 case 0
00277 typestr = 'ubit8';
00278 case 1
00279 typestr = 'short';
00280 case 2
00281 typestr = 'int';
00282 case 3
00283 typestr = 'float';
00284 otherwise
00285 err = ErrEval(FuncName,'Err_Cannot read this data type');
00286 return;
00287 end
00288
00289 %read .BRIK file
00290 BrikName = sprintf('%s.BRIK', BrikName);
00291
00292 fidBRIK = fopen(BrikName, 'rb', Opt.MachineFormat);
00293 if (fidBRIK < 0),
00294 err = ErrEval(FuncName,'Err_Could not open .BRIK file');
00295 return;
00296 end
00297
00298 numpix=Info.DATASET_DIMENSIONS(1)*Info.DATASET_DIMENSIONS(2);
00299 numslices=length(Opt.Slices);
00300 numframes=length(Opt.Frames);
00301
00302 if isallslices & isallframes
00303 V = fread(fidBRIK, (Info.DATASET_DIMENSIONS(1) .* Info.DATASET_DIMENSIONS(2) .* Info.DATASET_DIMENSIONS(3) .* Info.DATASET_RANK(2)) , typestr);
00304 else
00305 V=zeros(1,numpix*numslices*numframes);
00306 for k=1:numframes
00307 frame=Opt.Frames(k);
00308 if isallslices
00309 fseek(fidBRIK, numpix*Info.DATASET_DIMENSIONS(3)*(frame-1)*Info.TypeBytes, 'bof');
00310 istrt=1+numpix*numslices*(k-1);
00311 istp=istrt-1+numpix*numslices;
00312 V(istrt:istp)=fread(fidBRIK,numpix*numslices,typestr);
00313 else
00314 for j=1:numslices
00315 slice=Opt.Slices(j);
00316 fseek(fidBRIK, numpix*(slice-1+Info.DATASET_DIMENSIONS(3)*(frame-1))*Info.TypeBytes, 'bof');
00317 istrt=1+numpix*(j-1+numslices*(k-1));
00318 istp=istrt-1+numpix;
00319 V(istrt:istp)=fread(fidBRIK,numpix,typestr);
00320 end
00321 end
00322 end
00323 end
00324 fclose (fidBRIK);
00325
00326 %scale if required
00327 if (Opt.Scale),
00328 facs=Info.BRICK_FLOAT_FACS(Opt.Frames);
00329 iscl = find (facs); %find non zero scales
00330 NperBrik = Info.DATASET_DIMENSIONS(1) .* Info.DATASET_DIMENSIONS(2) .* numslices;
00331 if (~isempty(iscl)),
00332 for j=1:1:length(iscl),
00333 istrt = 1+ (iscl(j)-1).*NperBrik;
00334 istp = istrt + NperBrik - 1;
00335 V(istrt:istp) = V(istrt:istp) .* facs(iscl(j));
00336 end
00337 end
00338 end
00339
00340 %if required, reshape V
00341 if(~DoVect),
00342 V = reshape(V, Info.DATASET_DIMENSIONS(1), Info.DATASET_DIMENSIONS(2),...
00343 numslices, numframes);
00344 else
00345 V = reshape(V, Info.DATASET_DIMENSIONS(1).* Info.DATASET_DIMENSIONS(2).* numslices, numframes);
00346 end
00347
00348
00349 err = 0;
00350 return;
00351
00352 %test code
00353 clear all
00354 Opt.Slices=[2 8];
00355 Opt.Frames=[1 2];
00356
00357 BrikName='./r1_time+orig.BRIK';
00358 [err, V, Info, ErrMessage] = BrikLoad(BrikName);
00359 size(V)
00360 figure(1);
00361 k=0;
00362 for i=1:2
00363 for j=1:2
00364 k=k+1;
00365 subplot(2,2,k);
00366 imagesc(V(:,:,Opt.Slices(i),Opt.Frames(j)));
00367 if (exist('spectral')) colormap(spectral); end
00368 colorbar;
00369 end
00370 end
00371
00372 clear all
00373 Opt.Slices=[2 8];
00374 Opt.Frames=[1 2];
00375
00376
00377 BrikName='./r1_time+orig.BRIK';
00378 [err, V, Info, ErrMessage] = BrikLoad(BrikName, Opt);
00379 figure(2);
00380 k=0;
00381 for i=1:2
00382 for j=1:2
00383 k=k+1;
00384 subplot(2,2,k);
00385 imagesc(V(:,:,i,j)); colormap(spectral); colorbar;
00386 end
00387 end
00388
00389
