Home > vbmeg > external > bioelectromagnetism > avw_img_read.m

avw_img_read

PURPOSE ^

avw_img_read - read Analyze format data image (*.img)

SYNOPSIS ^

function [ avw, machine ] = avw_img_read(fileprefix,IMGorient,machine)

DESCRIPTION ^

 avw_img_read - read Analyze format data image (*.img)
 
 [ avw, machine ] = avw_img_read(filename)
 
 filename - filename of Analyze file
 
 (Modified by M. Sato)
 Do not change orientation of image data written in the file
 
 Returned values:
 
 avw.hdr - a struct with image data parameters.
 avw.img - a 3D matrix of image data (double precision).
 
 The 3D matrix of image is written in .img file as
  for z=1:NZ
     for y=1:NY
        for x=1:NX
           read(img(x,y,z))
        end
     end
  end
 
 Default ANALYZE coordinate system is Left-handed LAS
 
 X-Y plane is Transverse
 X-Z plane is Coronal
 Y-Z plane is Sagittal
 
 X axis runs from patient right (low X) to patient Left (high X)
 Y axis runs from posterior (low Y) to Anterior (high Y)
 Z axis runs from inferior (low Z) to Superior (high Z)
 
 See also: avw_hdr_read (called by this function),
           avw_view, avw_write, avw_img_write, avw_flip

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [ avw, machine ] = avw_img_read(fileprefix,IMGorient,machine)
0002 % avw_img_read - read Analyze format data image (*.img)
0003 %
0004 % [ avw, machine ] = avw_img_read(filename)
0005 %
0006 % filename - filename of Analyze file
0007 %
0008 % (Modified by M. Sato)
0009 % Do not change orientation of image data written in the file
0010 %
0011 % Returned values:
0012 %
0013 % avw.hdr - a struct with image data parameters.
0014 % avw.img - a 3D matrix of image data (double precision).
0015 %
0016 % The 3D matrix of image is written in .img file as
0017 %  for z=1:NZ
0018 %     for y=1:NY
0019 %        for x=1:NX
0020 %           read(img(x,y,z))
0021 %        end
0022 %     end
0023 %  end
0024 %
0025 % Default ANALYZE coordinate system is Left-handed LAS
0026 %
0027 % X-Y plane is Transverse
0028 % X-Z plane is Coronal
0029 % Y-Z plane is Sagittal
0030 %
0031 % X axis runs from patient right (low X) to patient Left (high X)
0032 % Y axis runs from posterior (low Y) to Anterior (high Y)
0033 % Z axis runs from inferior (low Z) to Superior (high Z)
0034 %
0035 % See also: avw_hdr_read (called by this function),
0036 %           avw_view, avw_write, avw_img_write, avw_flip
0037 %
0038 
0039 
0040 % $Revision: 1332 $ $Date:: 2011-02-24 15:57:20 +0900#$
0041 
0042 % Licence:  GNU GPL, no express or implied warranties
0043 % History:  05/2002, Darren.Weber@flinders.edu.au
0044 %                    The Analyze format is copyright
0045 %                    (c) Copyright, 1986-1995
0046 %                    Biomedical Imaging Resource, Mayo Foundation
0047 %           01/2003, Darren.Weber@flinders.edu.au
0048 %                    - adapted for matlab v5
0049 %                    - revised all orientation information and handling
0050 %                      after seeking further advice from AnalyzeDirect.com
0051 %           03/2003, Darren.Weber@flinders.edu.au
0052 %                    - adapted for -ve pixdim values (non standard Analyze)
0053 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0054 
0055 if ~exist('fileprefix','var'),
0056   msg = sprintf('...no input fileprefix - see help avw_img_read\n\n');
0057   error(msg);
0058 end
0059 if ~exist('IMGorient','var'), IMGorient = ''; end
0060 if ~exist('machine','var'), machine = 'ieee-le'; end
0061 
0062 if findstr('.hdr',fileprefix),
0063   fileprefix = strrep(fileprefix,'.hdr','');
0064 end
0065 if findstr('.img',fileprefix),
0066   fileprefix = strrep(fileprefix,'.img','');
0067 end
0068 
0069 % MAIN
0070 
0071 % Read the file header
0072 [ avw, machine ] = avw_hdr_read(fileprefix,machine);
0073 
0074 avw = read_image(avw,IMGorient,machine);
0075 
0076 return
0077 
0078 
0079 
0080 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0081 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0082 function [ avw ] = read_image(avw,IMGorient,machine)
0083 
0084 fid = fopen(sprintf('%s.img',avw.fileprefix),'r',machine);
0085 if fid < 0,
0086   msg = sprintf('...cannot open file %s.img\n\n',avw.fileprefix);
0087   error(msg);
0088 end
0089 
0090 % Modified by M. Sato
0091 %ver = '[$Revision: 1332 $]';
0092 %fprintf('\nAVW_IMG_READ [v%s]\n',ver(12:16));  tic;
0093 
0094 % short int bitpix;    /* Number of bits per pixel; 1, 8, 16, 32, or 64. */
0095 % short int datatype      /* Datatype for this image set */
0096 % /*Acceptable values for datatype are*/
0097 % _NONE             0
0098 % _UNKNOWN          0  /*Unknown data type*/
0099 % _BINARY           1  /*Binary             ( 1 bit per voxel)*/
0100 % _UNSIGNED_CHAR    2  /*Unsigned character ( 8 bits per voxel)*/
0101 % _SIGNED_SHORT     4  /*Signed short       (16 bits per voxel)*/
0102 % _SIGNED_INT       8  /*Signed integer     (32 bits per voxel)*/
0103 % _FLOAT           16  /*Floating point     (32 bits per voxel)*/
0104 % _COMPLEX         32  /*Complex (64 bits per voxel; 2 float)/*
0105 % _DOUBLE          64  /*Double precision   (64 bits per voxel)*/
0106 % _RGB            128  /*A Red-Green-Blue datatype*/
0107 % _ALL            255  /*Undocumented*/
0108 
0109 switch double(avw.hdr.dime.bitpix),
0110   case  1,   precision = 'bit1';
0111   case  8,   precision = 'uchar';
0112   case 16,   precision = 'int16';
0113   case 32,
0114     if     isequal(avw.hdr.dime.datatype, 8), precision = 'int32';
0115     else                                      precision = 'single';
0116     end
0117   case 64,   precision = 'double';
0118   otherwise,
0119     precision = 'uchar';
0120     fprintf('...precision undefined in header, using ''uchar''\n');
0121 end
0122 
0123 % read the whole .img file into matlab (faster)
0124 fprintf('...reading %s Analyze %s image format.\n',machine,precision);
0125 
0126 fseek(fid,0,'bof');
0127 
0128 % adjust for matlab version
0129 ver = version;
0130 ver = str2num(ver(1));
0131 if ver < 6,
0132   tmp = fread(fid,inf,sprintf('%s',precision));
0133 else,
0134   tmp = fread(fid,inf,sprintf('%s=>double',precision));
0135 end
0136 fclose(fid);
0137 
0138 % Update the global min and max values
0139 avw.hdr.dime.glmax = max(double(tmp));
0140 avw.hdr.dime.glmin = min(double(tmp));
0141 
0142 
0143 %---------------------------------------------------------------
0144 % Now partition the img data into xyz
0145 % dim[0]      Number of dimensions in database; usually 4.
0146 % dim[1]      Image X dimension;  number of pixels in an image row.
0147 % dim[2]      Image Y dimension;  number of pixel rows in slice.
0148 % dim[3]      Volume Z dimension; number of slices in a volume.
0149 % dim[4]      Time points; number of volumes in database.
0150 
0151 PixelDim = double(avw.hdr.dime.dim(2));
0152 RowDim   = double(avw.hdr.dime.dim(3));
0153 SliceDim = double(avw.hdr.dime.dim(4));
0154 
0155 PixelSz  = double(avw.hdr.dime.pixdim(2));
0156 RowSz    = double(avw.hdr.dime.pixdim(3));
0157 SliceSz  = double(avw.hdr.dime.pixdim(4));
0158 
0159 % ---- NON STANDARD ANALYZE...
0160 
0161 % Some Analyze files have been found to set -ve pixdim values, eg
0162 % the MNI template avg152T1_brain in the FSL etc/standard folder,
0163 % perhaps to indicate flipped orientation?  If so, this code below
0164 % will NOT handle the flip correctly!
0165 if PixelSz < 0,
0166   warning('X pixdim < 0 !!! resetting to abs(avw.hdr.dime.pixdim(2))');
0167   PixelSz = abs(PixelSz);
0168   avw.hdr.dime.pixdim(2) = single(PixelSz);
0169 end
0170 if RowSz < 0,
0171   warning('Y pixdim < 0 !!! resetting to abs(avw.hdr.dime.pixdim(3))');
0172   RowSz = abs(RowSz);
0173   avw.hdr.dime.pixdim(3) = single(RowSz);
0174 end
0175 if SliceSz < 0,
0176   warning('Z pixdim < 0 !!! resetting to abs(avw.hdr.dime.pixdim(4))');
0177   SliceSz = abs(SliceSz);
0178   avw.hdr.dime.pixdim(4) = single(SliceSz);
0179 end
0180 
0181 % ---- END OF NON STANDARD ANALYZE
0182 
0183 % --- Do not change orientation of image data
0184 %
0185     
0186 % --- Modified by M. Sato
0187 
0188 avw.img = reshape(tmp, [PixelDim, RowDim, SliceDim]);
0189     
0190 %    n = 1;
0191 %    x = 1:PixelDim;
0192 %    for z = 1:SliceDim,
0193 %      for y = 1:RowDim,
0194 %        % load Y row of X values into Z slice avw.img
0195 %        avw.img(x,y,z) = tmp(n:n+(PixelDim-1));
0196 %        n = n + PixelDim;
0197 %      end
0198 %    end
0199 
0200 % no need to rearrange avw.hdr.dime.dim or avw.hdr.dime.pixdim
0201     
0202 %
0203 %t=toc; fprintf('...done (%5.2f sec).\n\n',t);
0204 
0205 return
0206 
0207 % From ANALYZE documentation...
0208 %
0209 % The ANALYZE coordinate system has an origin in the lower left
0210 % corner. That is, with the subject lying supine, the coordinate
0211 % origin is on the right side of the body (x), at the back (y),
0212 % and at the feet (z). This means that:
0213 %
0214 % +X increases from right (R) to left (L)
0215 % +Y increases from the back (posterior,P) to the front (anterior, A)
0216 % +Z increases from the feet (inferior,I) to the head (superior, S)
0217 %
0218 % The LAS orientation is the radiological convention, where patient
0219 % left is on the image right.  The alternative neurological
0220 % convention is RAS (also Talairach convention).
0221 %
0222 % A major advantage of the Analzye origin convention is that the
0223 % coordinate origin of each orthogonal orientation (transverse,
0224 % coronal, and sagittal) lies in the lower left corner of the
0225 % slice as it is displayed.
0226 %
0227 % Orthogonal slices are numbered from one to the number of slices
0228 % in that orientation. For example, a volume (x, y, z) dimensioned
0229 % 128, 256, 48 has:
0230 %
0231 %   128 sagittal   slices numbered 1 through 128 (X)
0232 %   256 coronal    slices numbered 1 through 256 (Y)
0233 %    48 transverse slices numbered 1 through  48 (Z)
0234 %
0235 % Pixel coordinates are made with reference to the slice numbers from
0236 % which the pixels come. Thus, the first pixel in the volume is
0237 % referenced p(1,1,1) and not at p(0,0,0).
0238 %
0239 % Transverse slices are in the XY plane (also known as axial slices).
0240 % Sagittal slices are in the ZY plane.
0241 % Coronal slices are in the ZX plane.
0242 %
0243 %----------------------------------------------------------------------------

Generated on Mon 22-May-2023 06:53:56 by m2html © 2005