00001 function [err, ErrMessage, Info] = WriteBrik (M, Info, Opt)
00002 %
00003 % [err, ErrMessage, Info] = WriteBrik (M, Info, Opt)
00004 %
00005 %Purpose:
00006 %
00007 % Writes a data vector or matrix in AFNI's format
00008 % Also write data in AFNI's 1D file format
00009 %
00010 %Input Parameters:
00011 % M is the brick data (in vector or matrix format)
00012 % for 1D formats, M must be one or 2 dimensional.
00013 %
00014 % Info, the header structure (see BrikInfo)
00015 % if BYTEORDER_STRING is not specified, native is the default format
00016 % if BRICK_TYPES is not specified, short is the default.
00017 %
00018 % Info can be empty for 1D files but if you have one from
00019 % BrikLoad you should use it.
00020 %
00021 % Opt an options structure with the following fields, [default]
00022 % Most of the options are irrelevant for 1D formats.
00023 % .Scale: ([0]/1) scales values to a |maximum| of 32700.
00024 % This is useful for writing bricks as shorts.
00025 % .Prefix : filename prefix (mandatory field for 1D and BRIK formats)
00026 % .View : [+orig], +acpc, +tlrc
00027 % .verbose: ([0]/1)
00028 % .AppendHistory: (0/[1]) adds to the history field
00029 % .NoCheck: ([0]/1) flag that determines what checking should
00030 % be done to the Brick's header before writing.
00031 % 0: Full checking. This slows the function down
00032 % But you should use it whenever you're still
00033 % developing your code.
00034 % + The Header is passed through the function
00035 % CheckBrikHEAD
00036 % 1: + No Header Checking
00037 %
00038 % Regardless of .NoCheck, the following is done:
00039 % + The fields Info.BRICK_STATS and Info.BRICK_FLOAT_FACS
00040 % are cleared (and set if .Scale is used).
00041 % + The field Info.IDCODE_STRING is cleared in WriteBrikHEAD
00042 % + The field Info.IDCODE_DATE is set in WriteBrikHEAD
00043 %
00044 % .Slices: vector of slices, 1 based. Default is all slices.
00045 % .Slices can also be used to write a 1D dataset one chunk
00046 % at a time. For 1D files however, .Slices can only have 1 number.
00047 % When .Slices == 1 the file is open in write 'w' mode and M is
00048 % written into a new file.
00049 % When .Slices > 1 the file is open in append 'a' mode and M is
00050 % written into the end of the file.
00051 % .Frames: vector of frames, 1 based. Default is all frames.
00052 % If min(Slices)>1 or min(Frames)>1 then it assumes that the header
00053 % has already been written.
00054 % If Slices and Frames are not all slices and frames, then Scale
00055 % is set to 0 (no scaling).
00056 %
00057 %
00058 %Output Parameters:
00059 % err : 0 No Problem
00060 % : 1 Mucho Problems
00061 % ErrMessage : the error or warning message
00062 % Info : the header info of the written brick (or 1D file)
00063 %
00064 %Key Terms:
00065 %
00066 %More Info :
00067 % BrikInfo, BrikLoad, WriteBrikHEAD, HEAD_Rules, Info_1D
00068 %
00069 % need a FormHEAD function to create a minimal Info structure for a data vector
00070 %
00071 % version 2.0 (keep in sync. with WriteBrikHEAD)
00072 % In this version, Info.BRICK_STATS is cleared before writing.
00073 %
00074 % .Slices and .Frames options were added by Dr. Keith Worsley for FMRISTAT
00075 % http:
00076 % http:
00077 %
00078 % Author : Ziad Saad
00079 % Date : Fri Apr 6 16:03:02 PDT 2001, last modified Oct 01 04
00080 % LBC/NIMH/ National Institutes of Health, Bethesda Maryland
00081 %
00082 % Contact: ziad@nih.gov
00083
00084
00085 %Define the function name for easy referencing
00086 FuncName = 'WriteBrik';
00087 FUNCTION_VERSION = 'V2.0';
00088
00089 %Debug Flag
00090 DBG = 1;
00091
00092 %initailize return variables
00093 err = 1;
00094 ErrMessage = '';
00095
00096 if (nargin < 3),
00097 err = 1;
00098 ErrMessage = sprintf('Error %s: Need three input parameters', FuncName);
00099 errordlg(ErrMessage); return;
00100 end
00101
00102 %first check on Prefix
00103 if (~isfield(Opt, 'Prefix') | isempty (Opt.Prefix)),
00104 err = 1;
00105 ErrMessage = sprintf('Error %s: You must specify Opt.Prefix.', FuncName);
00106 errordlg(ErrMessage); return;
00107 end
00108 %set format if not present
00109 if (~isfield(Info,'FileFormat') | isempty(Info.FileFormat)),
00110 Info.FileFormat = 'BRIK';
00111 end
00112
00113 %is this a 1D file ?
00114 is1D = 0;
00115 if (isempty(Info) & size(M,4) == 1 & size(M,3) == 1),
00116 is1D = 1;
00117 elseif (strcmp(Info.FileFormat,'1D')),
00118 is1D = 1;
00119 if (size(M,4) ~= 1 | size(M,3) ~= 1),
00120 err = 1;
00121 ErrMessage = sprintf('Error %s: M is not one or two dimensional.', FuncName);
00122 errordlg(ErrMessage); return;
00123 end
00124 end
00125
00126
00127 if (is1D),
00128 if (isempty(Info)), Info = Info_1D(M); end
00129 Opt1D.OverWrite = 'n'; % default mode
00130 if (isfield(Opt,'Slices') & ~isempty(Opt.Slices)),
00131 if (length(Opt.Slices) ~= 1),
00132 err = 1;
00133 ErrMessage = sprintf('Error %s: .Slices can only have one value for 1D files', FuncName);
00134 errordlg(ErrMessage); return;
00135 end
00136 if (Opt.Slices > 1) Opt1D.OverWrite = 'a'; end
00137 end
00138 Opt1D.Space = 't';
00139 Opt1D.Fast = 'y';
00140 Opt1D.verbose = 0;
00141 [Name, Ext] = Remove1DExtension(Opt.Prefix);
00142 if (isempty(Ext)),
00143 if (~isempty(Info.Extension_1D)),
00144 Ext = sprintf('%s', Info.Extension_1D);
00145 else
00146 Ext = sprintf('.1D');
00147 end
00148 end
00149 FullName = sprintf('%s%s', Name, Ext);
00150 [err, UsedName] = wryte3(M, FullName, Opt1D);
00151 Info = Info_1D(M);
00152
00153 if (isempty(Ext) == 1),
00154 Ext = '.1D.dset';
00155 end
00156
00157 Info.Extension_1D = sprintf('%s', Ext);
00158 Info.RootName = sprintf('%s', Name);
00159 return;
00160 end
00161
00162 %check on options
00163 if (~isfield(Opt, 'verbose') | isempty (Opt.verbose)), Opt.verbose = 0; end
00164 if (~isfield(Opt, 'NoCheck') | isempty (Opt.NoCheck)), Opt.NoCheck = 0; end
00165
00166 if (Opt.verbose), fprintf(1,'%s verbose: Checking input data ...', FuncName); end
00167 if (~isfield(Opt, 'Scale') | isempty (Opt.Scale)), Opt.Scale = 0; end
00168 if (~isfield(Opt, 'View') | isempty(Opt.View)), Opt.View = 'orig'; end
00169 if (~isempty(findstr('orig', lower(Opt.View)))),
00170 Opt.Views = '+orig';
00171 elseif (~isempty(findstr('acpc', lower(Opt.View)))),
00172 Opt.Views = '+acpc';
00173 elseif (~isempty(findstr('tlrc', lower(Opt.View)))),
00174 Opt.Views = '+tlrc';
00175 else
00176 err = 1; ErrMessage = sprintf('Error %s: Bad value (%s) for Opt.View', FuncName, Opt.View); errordlg(ErrMessage); return;
00177 end
00178 if (~isfield(Opt, 'AppendHistory') | isempty (Opt.AppendHistory)), Opt.AppendHistory = 1; end
00179
00180 %form the flename based on the stuff in Opt.Prefix, just use the option
00181 Fname = sprintf('%s%s', Opt.Prefix, Opt.Views);
00182 FnameHEAD = sprintf('%s%s.HEAD', Opt.Prefix, Opt.Views);
00183 FnameBRIK = sprintf('%s%s.BRIK', Opt.Prefix, Opt.Views);
00184
00185 % This check is done later on before we write Slice 1 Frame 1 (see below)
00186 %if (exist(FnameHEAD) == 2 | exist(FnameBRIK) == 2),
00187 % err = 1; ErrMessage = sprintf('Error %s: Output data set %s exists.', FuncName, Fname); errordlg(ErrMessage); return;
00188 %end
00189
00190 %make sure there is no clash in input data
00191 %is M a 4D or 1D
00192 N = zeros(1,4);
00193 [N(1), N(2), N(3), N(4)] = size(M);
00194 nd = ndims(M);
00195
00196 % unsqueeze the array sizes (Keith's addition to fix writing time series, one slice at a time. Oct 01 04 )
00197 if nd==3 & length(Opt.Slices)==1 & length(Opt.Frames)>1
00198 N=[N(1) N(2) 1 N(3)];
00199 end
00200
00201 if (nd <= 2)
00202 isVect = 1;
00203 else
00204 isVect = 0;
00205 end
00206
00207 if (isfield(Info, 'DATASET_DIMENSIONS') & length(Info.DATASET_DIMENSIONS) < 3 & length(Info.DATASET_DIMENSIONS) > 0),
00208 err = 1; ErrMessage = sprintf('Error %s: If you specify DATASET_DIMENSIONS it must be a vector of three elements', FuncName); errordlg(ErrMessage); return;
00209 end
00210 if (isfield(Info, 'DATASET_RANK') & length(Info.DATASET_RANK) < 2),
00211 err = 1; ErrMessage = sprintf('Error %s: If you specify DATASET_RANK it must be a vector of two elements', FuncName); errordlg(ErrMessage); return;
00212 end
00213
00214 if ((~isfield(Info, 'DATASET_DIMENSIONS') | isempty(Info.DATASET_DIMENSIONS)) & isVect)
00215 err = 1; ErrMessage = sprintf('Error %s: If M is a vector, you must specify DATASET_DIMENSIONS in Info', FuncName); errordlg(ErrMessage); return;
00216 end
00217 if ((~isfield(Info, 'DATASET_RANK') | isempty(Info.DATASET_RANK)) & isVect)
00218 err = 1; ErrMessage = sprintf('Error %s: If M is a vector, you must specify DATASET_RANK in Info', FuncName); errordlg(ErrMessage); return;
00219 end
00220
00221 if (isfield(Info, 'DATASET_DIMENSIONS') & ~isempty(Info.DATASET_DIMENSIONS) & ~isVect)
00222 % if (N(1) ~= Info.DATASET_DIMENSIONS(1) | N(2) ~= Info.DATASET_DIMENSIONS(2) | N(3) ~= Info.DATASET_DIMENSIONS(3) | N(4) ~= Info.DATASET_RANK(2))
00223 if (N(1) ~= Info.DATASET_DIMENSIONS(1) | N(2) ~= Info.DATASET_DIMENSIONS(2) | N(3) > Info.DATASET_DIMENSIONS(3) | N(4) > Info.DATASET_RANK(2))
00224 err = 1; ErrMessage = sprintf('Error %s: Dimensions mismatch between dimensions of M and Info.DATASET_DIMENSIONS, Info.DATASET_RANK.', FuncName); errordlg(ErrMessage); return;
00225 end
00226 end
00227
00228 %OK, setup .DATASET_DIMENSIONS and .DATASET_RANK if needed
00229 if (~isfield(Info, 'DATASET_DIMENSIONS') | isempty(Info.DATASET_DIMENSIONS)),
00230 Info.DATASET_DIMENSIONS = N(1:3);
00231 end
00232 if (~isfield(Info, 'DATASET_RANK') | isempty(Info.DATASET_RANK)),
00233 Info.DATASET_RANK = [3 N(4)];
00234 end
00235
00236 %any Mandatory variables have now been set, check on the Header content
00237
00238 %Check out the values in Info
00239 if (~Opt.NoCheck),
00240 if (Opt.verbose), fprintf(1,'HEAD Info structure ...'); end
00241 [err, ErrMessage, Info] = CheckBrikHEAD(Info);
00242 if (err),
00243 ErrMessage = sprintf ('%s: Error in CheckBrikHEAD.\n(%s)', ErrMessage);
00244 return;
00245 end
00246 end
00247
00248 %reshape to a vector
00249 % if (~isVect | nd == 2),
00250 % M = reshape(M, prod(N), 1);
00251 %end
00252
00253 %Delete the Brick_Stats, let afni take care of them
00254 if (isfield(Info,'BRICK_STATS')), rmfield(Info,'BRICK_STATS'); end
00255 if (isfield(Info,'BRICK_FLOAT_FACS')), rmfield(Info,'BRICK_FLOAT_FACS'); end
00256
00257 %figure out the ouput format
00258 if (~isfield(Info, 'BRICK_TYPES') | isempty(Info.BRICK_TYPES)),
00259 B_type = 1; %short
00260 else
00261 B_type = Info.BRICK_TYPES;
00262 end
00263
00264 if (~isfield(Info, 'BYTEORDER_STRING') | isempty(Info.BYTEORDER_STRING)),
00265 % set the order based on the machine, used to be: ByteOrder = 'native'; prior to 17 Feb 04
00266 [c_c, mx_c, ed_c] = computer;
00267 if (ed_c(1) == 'L')
00268 ByteOrder = 'ieee-le'; %Little Endian
00269 Info.BYTEORDER_STRING = 'LSB_FIRST';
00270 elseif (ed_c(1) == 'B')
00271 ByteOrder = 'ieee-be'; %Little Endian
00272 Info.BYTEORDER_STRING = 'MSB_FIRST';
00273 else
00274 err = 1; ErrMessage = sprintf('Error %s: %s byte order is ambiguous.', FuncName, Info.BYTEORDER_STRING);
00275 return;
00276 end
00277 else
00278 if (~isempty(strmatch('MSB_FIRST', Info.BYTEORDER_STRING))),
00279 ByteOrder = 'ieee-be'; %Big Endian
00280 else
00281 if (~isempty(strmatch('LSB_FIRST', Info.BYTEORDER_STRING))),
00282 ByteOrder = 'ieee-le'; %Little Endian
00283 else
00284 err = 1; ErrMessage = sprintf('Error %s: %s byte order is ambiguous.', FuncName, Info.BYTEORDER_STRING);
00285 return;
00286 end
00287 end
00288 end
00289
00290 itype = unique(B_type);
00291 if (length(itype) > 1),
00292 err = 1; ErrMessage = sprintf('Error %s: Not set up to write sub-bricks of multiple sub-types', FuncName); errordlg(ErrMessage); return;
00293 end
00294 switch itype,
00295 case 0
00296 typestr = 'ubit8';
00297 case 1
00298 typestr = 'short';
00299 case 2
00300 typestr = 'int';
00301 case 3
00302 typestr = 'float';
00303 otherwise
00304 err = ErrEval(FuncName,'Err_Cannot write this data type');
00305 return;
00306 end
00307
00308 %figure out if scaling is required
00309 allslices=1:Info.DATASET_DIMENSIONS(3);
00310 if (~isfield(Opt,'Slices') | isempty(Opt.Slices)),
00311 Opt.Slices = allslices;
00312 end
00313 isallslices=all(ismember(allslices,Opt.Slices));
00314 allframes=1:Info.DATASET_RANK(2);
00315 if (~isfield(Opt,'Frames') | isempty(Opt.Frames)),
00316 Opt.Frames = allframes;
00317 end
00318 isallframes=all(ismember(allframes,Opt.Frames));
00319
00320 Info.BRICK_FLOAT_FACS = zeros(1,Info.DATASET_RANK(2));
00321 if (Opt.Scale & isallslices & isallframes),
00322 if (Opt.verbose), fprintf(1,'Scaling ...'); end
00323 NperBrik = Info.DATASET_DIMENSIONS(1) .* Info.DATASET_DIMENSIONS(2) .* Info.DATASET_DIMENSIONS(3);
00324 for (j=1:1:Info.DATASET_RANK(2)),
00325 istrt = 1+ (j-1).*NperBrik;
00326 istp = istrt + NperBrik - 1;
00327 [max1, imax1] = max(abs(M(istrt:istp)));
00328 Info.BRICK_STATS(2.*(j-1)+1)= min(M(istrt:istp));
00329 [max2, imax2] = max(M(istrt:istp));
00330 Info.BRICK_STATS(2.*j)= max2;
00331 Info.BRICK_FLOAT_FACS(j) = max1 ./ 32700; % a bit lower than 32000
00332 if (Info.BRICK_FLOAT_FACS(j) == 0)
00333 Info.BRICK_FLOAT_FACS(j) = 1;
00334 else
00335 M(istrt:istp) = M(istrt:istp) ./ Info.BRICK_FLOAT_FACS(j);
00336 end
00337 end
00338 end
00339
00340 numpix=Info.DATASET_DIMENSIONS(1)*Info.DATASET_DIMENSIONS(2);
00341 numslices=length(Opt.Slices);
00342 numframes=length(Opt.Frames);
00343
00344 %open file for writing based on the type specified in Info
00345 if Opt.Slices(1)==1 & Opt.Frames(1)==1
00346 if (exist(FnameHEAD) == 2 | exist(FnameBRIK) == 2),
00347 err = 1; ErrMessage = sprintf('Error %s: Output data set %s exists.', FuncName, Fname); errordlg(ErrMessage); return;
00348 end
00349 [fid, mess] = fopen (FnameBRIK, 'w', ByteOrder);
00350 if (fid < 0),
00351 err = 1; ErrMessage = sprintf('Error %s: Could not open %s for writing \n(%s)', FuncName, FnameBRIK, mess); errordlg(ErrMessage); return;
00352 end
00353 if ~(isallslices & isallframes)
00354 for frame=1:Info.DATASET_RANK(2)
00355 fwrite(fid,zeros(1,numpix*Info.DATASET_DIMENSIONS(3)),typestr);
00356 end
00357 end
00358 else
00359 [fid, mess] = fopen (FnameBRIK, 'r+', ByteOrder);
00360 if (fid < 0),
00361 err = 1; ErrMessage = sprintf('Error %s: Could not open %s for re-writing \n(%s)', FuncName, FnameBRIK, mess); errordlg(ErrMessage); return;
00362 end
00363 end
00364
00365 %write the file
00366 if (Opt.verbose), fprintf(1,'Writing %s to disk ...', FnameBRIK); end
00367 if isallslices & isallframes
00368 cnt = fwrite(fid, M, typestr);
00369 else
00370 cnt=0;
00371 for k=1:numframes
00372 frame=Opt.Frames(k);
00373 if isallslices
00374 fseek(fid, numpix*Info.DATASET_DIMENSIONS(3)*(frame-1)*Info.TypeBytes, 'bof');
00375 istrt=1+numpix*numslices*(k-1);
00376 istp=istrt-1+numpix*numslices;
00377 cnt=cnt+fwrite(fid,M(istrt:istp),typestr);
00378 else
00379 for j=1:numslices
00380 slice=Opt.Slices(j);
00381 fseek(fid, numpix*(slice-1+Info.DATASET_DIMENSIONS(3)*(frame-1))*Info.TypeBytes, 'bof');
00382 istrt=1+numpix*(j-1+numslices*(k-1));
00383 istp=istrt-1+numpix;
00384 cnt=cnt+fwrite(fid,M(istrt:istp),typestr);
00385 end
00386 end
00387 end
00388 end
00389
00390 if (cnt ~= prod(size(M))),
00391 err = 1; ErrMessage = sprintf('Error %s: Could not write all %d values to %s\n. Deleting %s ...', FuncName, FnameBRIK, prod(size(M)), FnameBRIK); errordlg(ErrMessage);
00392 fclose (fid);
00393 if (filexist(FnameBRIK) == 2),
00394 delete(FnameBRIK);
00395 end
00396 return;
00397 end
00398
00399 %close the file
00400 fclose (fid);
00401 [ST,I] = dbstack;
00402
00403 if Opt.Slices(1)==1 & Opt.Frames(1)==1
00404 %add the history note if needed
00405 if (Opt.AppendHistory),
00406 OptHist.AFNI_Format = 1;
00407 if isunix
00408 [tmp, OptHist.PerSig] = unix('whoami');
00409 %remove this annoying tset message (some bug ....)
00410 [err, snl, Nlines] = GetNextLine(OptHist.PerSig, 2);
00411 if (Nlines >= 2),
00412 [err, OptHist.PerSig] = GetNextLine(OptHist.PerSig,Nlines);
00413 end
00414 if (tmp),
00415 OptHist.PerSig = sprintf('DunnoWho');
00416 else
00417 OptHist.PerSig = zdeblank(OptHist.PerSig);
00418 end
00419 else
00420 OptHist.PerSig = sprintf('Not UNIX');
00421 end
00422 [err,S] = HistoryTrace (OptHist);
00423 if (~isfield(Info,'HISTORY_NOTE') |isempty(Info.HISTORY_NOTE)), Info.HISTORY_NOTE = ''; end
00424 Info.HISTORY_NOTE = sprintf('%s\\n%s', Info.HISTORY_NOTE, S);
00425 end
00426
00427 %call for the function to write the header
00428 if (Opt.verbose), fprintf(1,'Writing %s to disk ...', FnameHEAD); end
00429 [err, ErrMessage] = WriteBrikHEAD (FnameHEAD, Info);
00430 if (err),
00431 err = 1; ErrMessage = sprintf('Error %s: An error occurred in WriteBrikHEAD.\n%s', FuncName, ErrMessage); errordlg(ErrMessage); return;
00432 end
00433 if (Opt.verbose), fprintf(1,'Done.\n'); end
00434 end
00435
00436 err = 0;
00437 return;
00438
