avw_img_write - write Analyze image files (*.img)
avw_img_write(avw, fileprefix, [IMGorient], [machine])
avw.img - a 3D matrix of image data (double precision).
avw.hdr - a struct with image data parameters. If
not empty, this function calls avw_hdr_write.
fileprefix - a string, the filename without the .img
extension. If empty, may use avw.fileprefix
IMGorient - optional int, force writing of specified
orientation, with values:
[], if empty, will use avw.hdr.hist.orient field
0, transverse/axial unflipped (default, radiological)
1, coronal unflipped
2, sagittal unflipped
3, transverse/axial flipped, left to right
4, coronal flipped, anterior to posterior
5, sagittal flipped, superior to inferior
This function will set avw.hdr.hist.orient and write the
image data in a corresponding order. This function is
in alpha development, so it has not been exhaustively
tested (07/2003). See avw_img_read for more information
and documentation on the orientation option.
Orientations 3-5 are NOT recommended! They are part
of the Analyze format, but only used in Analyze
for faster raster graphics during movies.
machine - a string, see machineformat in fread for details.
The default here is 'ieee-le'.
Tip: to change the data type, set avw.hdr.dime.datatype to:
1 Binary ( 1 bit per voxel)
2 Unsigned character ( 8 bits per voxel)
4 Signed short ( 16 bits per voxel)
8 Signed integer ( 32 bits per voxel)
16 Floating point ( 32 bits per voxel)
32 Complex, 2 floats ( 64 bits per voxel), flipping not supported
64 Double precision ( 64 bits per voxel)
128 Red-Green-Blue (128 bits per voxel), not supported
See also: avw_write, avw_hdr_write,
avw_read, avw_hdr_read, avw_img_read, avw_view0001 function avw_img_write(avw, fileprefix, IMGorient, machine) 0002 0003 % avw_img_write - write Analyze image files (*.img) 0004 % 0005 % avw_img_write(avw, fileprefix, [IMGorient], [machine]) 0006 % 0007 % avw.img - a 3D matrix of image data (double precision). 0008 % avw.hdr - a struct with image data parameters. If 0009 % not empty, this function calls avw_hdr_write. 0010 % 0011 % fileprefix - a string, the filename without the .img 0012 % extension. If empty, may use avw.fileprefix 0013 % 0014 % IMGorient - optional int, force writing of specified 0015 % orientation, with values: 0016 % 0017 % [], if empty, will use avw.hdr.hist.orient field 0018 % 0, transverse/axial unflipped (default, radiological) 0019 % 1, coronal unflipped 0020 % 2, sagittal unflipped 0021 % 3, transverse/axial flipped, left to right 0022 % 4, coronal flipped, anterior to posterior 0023 % 5, sagittal flipped, superior to inferior 0024 % 0025 % This function will set avw.hdr.hist.orient and write the 0026 % image data in a corresponding order. This function is 0027 % in alpha development, so it has not been exhaustively 0028 % tested (07/2003). See avw_img_read for more information 0029 % and documentation on the orientation option. 0030 % Orientations 3-5 are NOT recommended! They are part 0031 % of the Analyze format, but only used in Analyze 0032 % for faster raster graphics during movies. 0033 % 0034 % machine - a string, see machineformat in fread for details. 0035 % The default here is 'ieee-le'. 0036 % 0037 % Tip: to change the data type, set avw.hdr.dime.datatype to: 0038 % 0039 % 1 Binary ( 1 bit per voxel) 0040 % 2 Unsigned character ( 8 bits per voxel) 0041 % 4 Signed short ( 16 bits per voxel) 0042 % 8 Signed integer ( 32 bits per voxel) 0043 % 16 Floating point ( 32 bits per voxel) 0044 % 32 Complex, 2 floats ( 64 bits per voxel), flipping not supported 0045 % 64 Double precision ( 64 bits per voxel) 0046 % 128 Red-Green-Blue (128 bits per voxel), not supported 0047 % 0048 % See also: avw_write, avw_hdr_write, 0049 % avw_read, avw_hdr_read, avw_img_read, avw_view 0050 % 0051 0052 % $Revision: 1.1 $ $Date: 2004/11/12 01:30:25 $ 0053 0054 % Licence: GNU GPL, no express or implied warranties 0055 % History: 05/2002, Darren.Weber@flinders.edu.au 0056 % The Analyze format is copyright 0057 % (c) Copyright, 1986-1995 0058 % Biomedical Imaging Resource, Mayo Foundation 0059 % 07/2004, chodkowski@kennedykrieger.org, added ability to 0060 % write complex AVW .img files. added error if complex 0061 % data is to be flipped. currently there is no logic for flipping 0062 % complex data. also force invalid IMGorient to 0, ie, no flipping. 0063 % forcing an invalid IMGorient to 0 allows me to remove duplicate code 0064 % in the IMGorient case/otherwise logic. i also pulled the fwrite of 0065 % non-flipped data, IMGorient == 0, out of any looping mechanism. looping 0066 % is not necessary as the data is already in its correct orientation. 0067 % removing fwrite from the loops should be faster but, more importantly, 0068 % it allows the writing of complex data which contains two matrix 0069 % elements per one pixel. moved the fclose statement from write_image 0070 % to avw_img_write, the same function that calls fopen. 0071 % 0072 % write complex data example: 0073 % % assume rr contains n-dimensional real values 0074 % % assume ii contains n-dimensional imaginary values 0075 % % where the magnitude image can be computed as follows: 0076 % 0077 % magn = sqrt( rr.^2 + ii.^2 ); 0078 % 0079 % cc = zeros( prod( size( rr ))*2, 1 ); % "*2" bc 1pix = [real, imag] 0080 % cc(1:2:end) = reshape( rr, prod(size(rr)), 1 ); 0081 % cc(2:2:end) = reshape( ii, prod(size(ii)), 1 ); 0082 % 0083 % avw = avw_hdr_make; 0084 % 0085 % % numDims, xdim, ydim, zdim, ..., padded with 8 0s 0086 % tmp = [ ndims(rr) size(rr) zeros(1,8) ]; 0087 % avw.hdr.dime.dim = tmp(1:8); % keep only the first 8 values 0088 % 0089 % avw.hdr.dime.datatype = 32; % complex 0090 % avw.hdr.dime.bitpix = 64; % 4+4 bytes/pixel 0091 % 0092 % avw.img = cc; 0093 % avw_img_write( avw, 'foo' ); 0094 % 0095 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0096 0097 0098 %------------------------------------------------------------------------ 0099 % Check inputs 0100 0101 if ~exist('avw','var'), 0102 doc avw_img_write; 0103 msg = sprintf('\n...no input avw.\n'); 0104 error(msg); 0105 elseif isempty(avw), 0106 msg = sprintf('\n...empty input avw.\n'); 0107 error(msg); 0108 elseif ~isfield(avw,'img'), 0109 msg = sprintf('\n...empty input avw.img\n'); 0110 error(msg); 0111 end 0112 0113 if ~exist('fileprefix','var'), 0114 if isfield(avw,'fileprefix'), 0115 if ~isempty(avw.fileprefix), 0116 fileprefix = avw.fileprefix; 0117 else 0118 fileprefix = ''; 0119 end 0120 end 0121 end 0122 if isempty(fileprefix), 0123 doc avw_img_write; 0124 fprintf('\n...no input fileprefix - see help avw_img_write\n'); 0125 return; 0126 end 0127 0128 if ~exist('IMGorient','var'), IMGorient = ''; end 0129 if ~exist('machine','var'), machine = 'ieee-le'; end 0130 0131 if findstr('.hdr',fileprefix), 0132 % fprintf('AVW_IMG_WRITE: Removing .hdr extension from ''%s''\n',fileprefix); 0133 fileprefix = strrep(fileprefix,'.hdr',''); 0134 end 0135 if findstr('.img',fileprefix), 0136 % fprintf('AVW_IMG_WRITE: Removing .img extension from ''%s''\n',fileprefix); 0137 fileprefix = strrep(fileprefix,'.img',''); 0138 end 0139 0140 0141 0142 %------------------------------------------------------------------------ 0143 % MAIN 0144 0145 version = '[$Revision: 1.1 $]'; 0146 fprintf('\nAVW_IMG_WRITE [v%s]\n',version(12:16)); tic; 0147 0148 fid = fopen(sprintf('%s.img',fileprefix),'w',machine); 0149 if fid < 0, 0150 msg = sprintf('Cannot open file %s.img\n',fileprefix); 0151 error(msg); 0152 else 0153 avw = write_image(fid,avw,fileprefix,IMGorient,machine); 0154 end 0155 fclose(fid); 0156 0157 t=toc; fprintf('...done (%5.2f sec).\n\n',t); 0158 0159 % MUST write header after the image, to ensure any 0160 % orientation changes during image write are saved 0161 % in the header 0162 avw_hdr_write(avw,fileprefix,machine); 0163 0164 return 0165 0166 0167 0168 0169 %----------------------------------------------------------------------------------- 0170 0171 function avw = write_image(fid,avw,fileprefix,IMGorient,machine) 0172 0173 % short int bitpix; /* Number of bits per pixel; 1, 8, 16, 32, or 64. */ 0174 % short int datatype /* Datatype for this image set */ 0175 % /*Acceptable values for datatype are*/ 0176 % #define DT_NONE 0 0177 % #define DT_UNKNOWN 0 /*Unknown data type*/ 0178 % #define DT_BINARY 1 /*Binary ( 1 bit per voxel)*/ 0179 % #define DT_UNSIGNED_CHAR 2 /*Unsigned character ( 8 bits per voxel)*/ 0180 % #define DT_SIGNED_SHORT 4 /*Signed short (16 bits per voxel)*/ 0181 % #define DT_SIGNED_INT 8 /*Signed integer (32 bits per voxel)*/ 0182 % #define DT_FLOAT 16 /*Floating point (32 bits per voxel)*/ 0183 % #define DT_COMPLEX 32 /*Complex,2 floats (64 bits per voxel)/* 0184 % #define DT_DOUBLE 64 /*Double precision (64 bits per voxel)*/ 0185 % #define DT_RGB 128 /*A Red-Green-Blue datatype*/ 0186 % #define DT_ALL 255 /*Undocumented*/ 0187 0188 switch double(avw.hdr.dime.datatype), 0189 case 1, 0190 avw.hdr.dime.bitpix = int16( 1); precision = 'bit1'; 0191 case 2, 0192 avw.hdr.dime.bitpix = int16( 8); precision = 'uchar'; 0193 case 4, 0194 avw.hdr.dime.bitpix = int16(16); precision = 'int16'; 0195 case 8, 0196 avw.hdr.dime.bitpix = int16(32); precision = 'int32'; 0197 case 16, 0198 avw.hdr.dime.bitpix = int16(32); precision = 'single'; 0199 case 32, 0200 avw.hdr.dime.bitpix = int16(64); precision = 'float'; 0201 case 64, 0202 avw.hdr.dime.bitpix = int16(64); precision = 'double'; 0203 case 128, 0204 error('...RGB datatype not yet supported.\n'); 0205 otherwise 0206 fprintf('...unknown datatype, using type 16 (32 bit floats).\n'); 0207 avw.hdr.dime.datatype = int16(16); 0208 avw.hdr.dime.bitpix = int16(32); precision = 'single'; 0209 end 0210 0211 0212 % write the .img file, depending on the .img orientation 0213 fprintf('...writing %s precision Analyze image (%s).\n',precision,machine); 0214 0215 fseek(fid,0,'bof'); 0216 0217 % The standard image orientation is axial unflipped 0218 if isempty(avw.hdr.hist.orient), 0219 msg = [ '...WARNING: avw.hdr.hist.orient ~= 0.\n',... 0220 ' This function assumes the input avw.img is\n',... 0221 ' in axial unflipped orientation in memory. This is\n',... 0222 ' created by the avw_img_read function, which converts\n',... 0223 ' any input file image to axial unflipped in memory.\n']; 0224 fprintf(msg) 0225 end 0226 0227 if isempty(IMGorient), 0228 fprintf('...no IMGorient specified, using avw.hdr.hist.orient value.\n'); 0229 IMGorient = double(avw.hdr.hist.orient); 0230 end 0231 0232 if ~isfinite(IMGorient), 0233 fprintf('...WARNING: IMGorient is not finite!\n'); 0234 fprintf('...unknown orientation specified, assuming default axial unflipped\n'); 0235 fprintf('...using case 0 \n' ); 0236 IMGorient = 0; 0237 end 0238 0239 maxCase = 5; 0240 if (IMGorient > maxCase ), 0241 fprintf('...WARNING: IMGorient is greater than %d!\n', maxCase); 0242 fprintf('...unknown orientation specified, assuming default axial unflipped\n'); 0243 fprintf('...using case 0 \n' ); 0244 IMGorient = 0; 0245 end 0246 0247 % if datatype is complex (32), do not allow flipping bc IMGorient == 0 is 0248 % the only fwrite logic that will handle 2 elements/pixel. i modified the 0249 % non-flipping fwrite logic to write out then entire avw.img in a single 0250 % block, ie, fwrite does not care about the volume dimensions. this allows 0251 % one to write out complex data where one pixel is represented by 2 elements 0252 % in the avw.img matrix. note that AVW complex data means: 0253 % [ [r1,i1], [r2,i2], ... ] 0254 0255 % issue an error if complex data needs to be flipped 0256 if (( avw.hdr.dime.datatype == 32 ) && ( IMGorient ~= 0 )) 0257 msg = [ '...ERROR: avw.hdr.dime.datatype = 32 (complex) and IMGorient, ', ... 0258 'the orientation of the volume, is not 0. A non-zero ', ... 0259 'IMGorient requires flipping the data. Flipping is not ', ... 0260 'implemented for complex data. Flip your data before ', ... 0261 'calling this function' ]; 0262 msg = sprintf( '%s (%s).', msg, mfilename ); 0263 error( msg ); 0264 end; 0265 0266 0267 switch IMGorient, 0268 0269 case 0, % transverse/axial unflipped 0270 0271 % For the 'transverse unflipped' type, the voxels are stored with 0272 % Pixels in 'x' axis (varies fastest) - from patient right to left 0273 % Rows in 'y' axis - from patient posterior to anterior 0274 % Slices in 'z' axis - from patient inferior to superior 0275 0276 fprintf('...writing axial unflipped\n'); 0277 0278 avw.hdr.hist.orient = uint8(0); 0279 0280 SliceDim = double(avw.hdr.dime.dim(4)); % z 0281 RowDim = double(avw.hdr.dime.dim(3)); % y 0282 PixelDim = double(avw.hdr.dime.dim(2)); % x 0283 SliceSz = double(avw.hdr.dime.pixdim(4)); 0284 RowSz = double(avw.hdr.dime.pixdim(3)); 0285 PixelSz = double(avw.hdr.dime.pixdim(2)); 0286 0287 newWay = true; 0288 if ( newWay ) 0289 % since there is no flipping, write out the entire volume at once 0290 fwrite( fid, avw.img, precision ); 0291 else 0292 x = 1:PixelDim; 0293 for z = 1:SliceDim, 0294 for y = 1:RowDim, 0295 fwrite(fid,avw.img(x,y,z),precision); 0296 end 0297 end 0298 end; 0299 0300 case 1, % coronal unflipped 0301 0302 % For the 'coronal unflipped' type, the voxels are stored with 0303 % Pixels in 'x' axis (varies fastest) - from patient right to left 0304 % Rows in 'z' axis - from patient inferior to superior 0305 % Slices in 'y' axis - from patient posterior to anterior 0306 0307 fprintf('...writing coronal unflipped\n'); 0308 0309 avw.hdr.hist.orient = uint8(1); 0310 0311 SliceDim = double(avw.hdr.dime.dim(3)); % y 0312 RowDim = double(avw.hdr.dime.dim(4)); % z 0313 PixelDim = double(avw.hdr.dime.dim(2)); % x 0314 SliceSz = double(avw.hdr.dime.pixdim(3)); 0315 RowSz = double(avw.hdr.dime.pixdim(4)); 0316 PixelSz = double(avw.hdr.dime.pixdim(2)); 0317 0318 x = 1:PixelDim; 0319 for y = 1:SliceDim, 0320 for z = 1:RowDim, 0321 fwrite(fid,avw.img(x,y,z),precision); 0322 end 0323 end 0324 0325 case 2, % sagittal unflipped 0326 0327 % For the 'sagittal unflipped' type, the voxels are stored with 0328 % Pixels in 'y' axis (varies fastest) - from patient posterior to anterior 0329 % Rows in 'z' axis - from patient inferior to superior 0330 % Slices in 'x' axis - from patient right to left 0331 0332 fprintf('...writing sagittal unflipped\n'); 0333 0334 avw.hdr.hist.orient = uint8(2); 0335 0336 SliceDim = double(avw.hdr.dime.dim(2)); % x 0337 RowDim = double(avw.hdr.dime.dim(4)); % z 0338 PixelDim = double(avw.hdr.dime.dim(3)); % y 0339 SliceSz = double(avw.hdr.dime.pixdim(2)); 0340 RowSz = double(avw.hdr.dime.pixdim(4)); 0341 PixelSz = double(avw.hdr.dime.pixdim(3)); 0342 0343 y = 1:PixelDim; 0344 for x = 1:SliceDim, 0345 for z = 1:RowDim, 0346 fwrite(fid,avw.img(x,y,z),precision); 0347 end 0348 end 0349 0350 case 3, % transverse/axial flipped 0351 0352 % For the 'transverse flipped' type, the voxels are stored with 0353 % Pixels in 'x' axis (varies fastest) - from patient right to left 0354 % Rows in 'y' axis - from patient anterior to posterior* 0355 % Slices in 'z' axis - from patient inferior to superior 0356 0357 fprintf('...writing axial flipped (+Y from Anterior to Posterior)\n'); 0358 0359 avw.hdr.hist.orient = uint8(3); 0360 0361 SliceDim = double(avw.hdr.dime.dim(4)); % z 0362 RowDim = double(avw.hdr.dime.dim(3)); % y 0363 PixelDim = double(avw.hdr.dime.dim(2)); % x 0364 SliceSz = double(avw.hdr.dime.pixdim(4)); 0365 RowSz = double(avw.hdr.dime.pixdim(3)); 0366 PixelSz = double(avw.hdr.dime.pixdim(2)); 0367 0368 x = 1:PixelDim; 0369 for z = 1:SliceDim, 0370 for y = RowDim:-1:1, % flipped in Y 0371 fwrite(fid,avw.img(x,y,z),precision); 0372 end 0373 end 0374 0375 case 4, % coronal flipped 0376 0377 % For the 'coronal flipped' type, the voxels are stored with 0378 % Pixels in 'x' axis (varies fastest) - from patient right to left 0379 % Rows in 'z' axis - from patient inferior to superior 0380 % Slices in 'y' axis - from patient anterior to posterior 0381 0382 fprintf('...writing coronal flipped (+Z from Superior to Inferior)\n'); 0383 0384 avw.hdr.hist.orient = uint8(4); 0385 0386 SliceDim = double(avw.hdr.dime.dim(3)); % y 0387 RowDim = double(avw.hdr.dime.dim(4)); % z 0388 PixelDim = double(avw.hdr.dime.dim(2)); % x 0389 SliceSz = double(avw.hdr.dime.pixdim(3)); 0390 RowSz = double(avw.hdr.dime.pixdim(4)); 0391 PixelSz = double(avw.hdr.dime.pixdim(2)); 0392 0393 x = 1:PixelDim; 0394 for y = 1:SliceDim, 0395 for z = RowDim:-1:1, 0396 fwrite(fid,avw.img(x,y,z),precision); 0397 end 0398 end 0399 0400 case 5, % sagittal flipped 0401 0402 % For the 'sagittal flipped' type, the voxels are stored with 0403 % Pixels in 'y' axis (varies fastest) - from patient posterior to anterior 0404 % Rows in 'z' axis - from patient superior to inferior 0405 % Slices in 'x' axis - from patient right to left 0406 0407 fprintf('...writing sagittal flipped (+Z from Superior to Inferior)\n'); 0408 0409 avw.hdr.hist.orient = uint8(5); 0410 0411 SliceDim = double(avw.hdr.dime.dim(2)); % x 0412 RowDim = double(avw.hdr.dime.dim(4)); % z 0413 PixelDim = double(avw.hdr.dime.dim(3)); % y 0414 SliceSz = double(avw.hdr.dime.pixdim(2)); 0415 RowSz = double(avw.hdr.dime.pixdim(4)); 0416 PixelSz = double(avw.hdr.dime.pixdim(3)); 0417 0418 y = 1:PixelDim; 0419 for x = 1:SliceDim, 0420 for z = RowDim:-1:1, % superior to inferior 0421 fwrite(fid,avw.img(x,y,z),precision); 0422 end 0423 end 0424 0425 0426 %% this fall-thru case should never happen since an unknown IMGorient 0427 %% is detected before this case statement and forced to "0". no need 0428 %% to have duplicate logic as this block of code is the same as case 0. 0429 0430 %otherwise, % transverse/axial unflipped 0431 % 0432 % % For the 'transverse unflipped' type, the voxels are stored with 0433 % % Pixels in 'x' axis (varies fastest) - from patient right to left 0434 % % Rows in 'y' axis - from patient posterior to anterior 0435 % % Slices in 'z' axis - from patient inferior to superior 0436 % 0437 % fprintf('...unknown orientation specified, assuming default axial unflipped\n'); 0438 % 0439 % avw.hdr.hist.orient = uint8(0); 0440 % 0441 % SliceDim = double(avw.hdr.dime.dim(4)); % z 0442 % RowDim = double(avw.hdr.dime.dim(3)); % y 0443 % PixelDim = double(avw.hdr.dime.dim(2)); % x 0444 % SliceSz = double(avw.hdr.dime.pixdim(4)); 0445 % RowSz = double(avw.hdr.dime.pixdim(3)); 0446 % PixelSz = double(avw.hdr.dime.pixdim(2)); 0447 % 0448 % x = 1:PixelDim; 0449 % for z = 1:SliceDim, 0450 % for y = 1:RowDim, 0451 % fwrite(fid,avw.img(x,y,z),precision); 0452 % end 0453 % end 0454 0455 end 0456 0457 % Update the header 0458 avw.hdr.dime.dim(2:4) = int16([PixelDim,RowDim,SliceDim]); 0459 avw.hdr.dime.pixdim(2:4) = single([PixelSz,RowSz,SliceSz]); 0460 0461 return