avw_img_read.m

mri_toolbox是一个工具用来MRI. 来自于SourceForge, 我上传这个软件,希望能结识对医疗软件感兴趣的兄弟.

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

% avw_img_read - read Analyze format data image (*.img)
% 
% [ avw, machine ] = avw_img_read(fileprefix, [orient], [machine])
% 
% fileprefix - a string, the filename without the .img extension
% 
% orient - read a specified orientation, integer values:
% 
%          '', use header history orient field
%          0,  transverse unflipped (LAS*)
%          1,  coronal unflipped (LA*S)
%          2,  sagittal unflipped (L*AS)
%          3,  transverse flipped (LPS*)
%          4,  coronal flipped (LA*I)
%          5,  sagittal flipped (L*AI)
% 
% where * follows the slice dimension and letters indicate +XYZ
% orientations (L left, R right, A anterior, P posterior,
% I inferior, & S superior).
% 
% Some files may contain data in the 3-5 orientations, but this
% is unlikely. For more information about orientation, see the
% documentation at the end of this .m file.  See also the 
% AVW_FLIP function for orthogonal reorientation.
% 
% machine - a string, see machineformat in fread for details.
%           The default here is 'ieee-le' but the routine
%           will automatically switch between little and big
%           endian to read any such Analyze header.  It
%           reports the appropriate machine format and can
%           return the machine value.
% 
% Returned values:
% 
% avw.hdr - a struct with image data parameters.
% avw.img - a 3D matrix of image data (double precision).
% 
% The returned 3D matrix will correspond with the 
% default ANALYZE coordinate system, which 
% is Left-handed:
% 
% 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
% 


% $Revision: 1.9 $ $Date: 2004/02/07 01:41:51 $

% Licence:  GNU GPL, no express or implied warranties
% History:  05/2002, [e-mail]
%                    The Analyze format is copyright 
%                    (c) Copyright, 1986-1995
%                    Biomedical Imaging Resource, Mayo Foundation
%           01/2003, [e-mail]
%                    - adapted for matlab v5
%                    - revised all orientation information and handling 
%                      after seeking further advice from AnalyzeDirect.com
%           03/2003, [e-mail]
%                    - adapted for -ve pixdim values (non standard Analyze)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('fileprefix','var'),
  msg = sprintf('...no input fileprefix - see help avw_img_read\n\n');
  error(msg);
end
if ~exist('IMGorient','var'), IMGorient = ''; end
if ~exist('machine','var'), machine = 'ieee-le'; end

if findstr('.hdr',fileprefix),
  fileprefix = strrep(fileprefix,'.hdr','');
end
if findstr('.img',fileprefix),
  fileprefix = strrep(fileprefix,'.img','');
end

% MAIN

% Read the file header
[ avw, machine ] = avw_hdr_read(fileprefix,machine);

avw = read_image(avw,IMGorient,machine);

return



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ avw ] = read_image(avw,IMGorient,machine)

fid = fopen(sprintf('%s.img',avw.fileprefix),'r',machine);
if fid < 0,
  msg = sprintf('...cannot open file %s.img\n\n',avw.fileprefix);
  error(msg);
end

ver = '[$Revision: 1.9 $]';
fprintf('\nAVW_IMG_READ [v%s]\n',ver(12:16));  tic;

% short int bitpix;    /* Number of bits per pixel; 1, 8, 16, 32, or 64. */ 
% short int datatype      /* Datatype for this image set */ 
% /*Acceptable values for datatype are*/ 
% #define DT_NONE             0
% #define DT_UNKNOWN          0    /*Unknown data type*/ 
% #define DT_BINARY           1    /*Binary             ( 1 bit per voxel)*/ 
% #define DT_UNSIGNED_CHAR    2    /*Unsigned character ( 8 bits per voxel)*/ 
% #define DT_SIGNED_SHORT     4    /*Signed short       (16 bits per voxel)*/ 
% #define DT_SIGNED_INT       8    /*Signed integer     (32 bits per voxel)*/ 
% #define DT_FLOAT           16    /*Floating point     (32 bits per voxel)*/ 
% #define DT_COMPLEX         32    /*Complex (64 bits per voxel; 2 floating point numbers)/* 
% #define DT_DOUBLE          64    /*Double precision   (64 bits per voxel)*/ 
% #define DT_RGB            128    /*A Red-Green-Blue datatype*/
% #define DT_ALL            255    /*Undocumented*/

switch double(avw.hdr.dime.bitpix),
  case  1,   precision = 'bit1';
  case  8,   precision = 'uchar';
  case 16,   precision = 'int16';
  case 32,
    if     isequal(avw.hdr.dime.datatype, 8), precision = 'int32';
    else                                      precision = 'single';
    end
  case 64,   precision = 'double';
  otherwise,
    precision = 'uchar';
    fprintf('...precision undefined in header, using ''uchar''\n');
end

% read the whole .img file into matlab (faster)
fprintf('...reading %s Analyze %s image format.\n',machine,precision);
fseek(fid,0,'bof');
% adjust for matlab version
ver = version;
ver = str2num(ver(1));
if ver < 6,
  tmp = fread(fid,inf,sprintf('%s',precision));
else,
  tmp = fread(fid,inf,sprintf('%s=>double',precision));
end
fclose(fid);

% Update the global min and max values
avw.hdr.dime.glmax = max(double(tmp));
avw.hdr.dime.glmin = min(double(tmp));


%---------------------------------------------------------------
% Now partition the img data into xyz

% --- first figure out the size of the image

% short int dim[ ];      /* Array of the image dimensions */ 
%
% dim[0]      Number of dimensions in database; usually 4. 
% dim[1]      Image X dimension;  number of pixels in an image row. 
% dim[2]      Image Y dimension;  number of pixel rows in slice. 
% dim[3]      Volume Z dimension; number of slices in a volume. 
% dim[4]      Time points; number of volumes in database.

PixelDim = double(avw.hdr.dime.dim(2));
RowDim   = double(avw.hdr.dime.dim(3));
SliceDim = double(avw.hdr.dime.dim(4));

PixelSz  = double(avw.hdr.dime.pixdim(2));
RowSz    = double(avw.hdr.dime.pixdim(3));
SliceSz  = double(avw.hdr.dime.pixdim(4));





% ---- NON STANDARD ANALYZE...

% Some Analyze files have been found to set -ve pixdim values, eg
% the MNI template avg152T1_brain in the FSL etc/standard folder,
% perhaps to indicate flipped orientation?  If so, this code below
% will NOT handle the flip correctly!
if PixelSz < 0,
  warning('X pixdim < 0 !!! resetting to abs(avw.hdr.dime.pixdim(2))');
  PixelSz = abs(PixelSz);
  avw.hdr.dime.pixdim(2) = single(PixelSz);
end
if RowSz < 0,
  warning('Y pixdim < 0 !!! resetting to abs(avw.hdr.dime.pixdim(3))');
  RowSz = abs(RowSz);
  avw.hdr.dime.pixdim(3) = single(RowSz);
end
if SliceSz < 0,
  warning('Z pixdim < 0 !!! resetting to abs(avw.hdr.dime.pixdim(4))');
  SliceSz = abs(SliceSz);
  avw.hdr.dime.pixdim(4) = single(SliceSz);
end

% ---- END OF NON STANDARD ANALYZE





% --- check the orientation specification and arrange img accordingly
if ~isempty(IMGorient),
  if ischar(IMGorient),
    avw.hdr.hist.orient = uint8(str2num(IMGorient));
  else
    avw.hdr.hist.orient = uint8(IMGorient);
  end
end,

if isempty(avw.hdr.hist.orient),
  msg = [ '...unspecified avw.hdr.hist.orient, using default 0\n',...
      '   (check image and try explicit IMGorient option).\n'];
  fprintf(msg);
  avw.hdr.hist.orient = uint8(0);
end

switch double(avw.hdr.hist.orient),
  
  case 0, % transverse unflipped
    
    % orient = 0:  The primary orientation of the data on disk is in the
    % transverse plane relative to the object scanned.  Most commonly, the fastest
    % moving index through the voxels that are part of this transverse image would
    % span the right-left extent of the structure imaged, with the next fastest
    % moving index spanning the posterior-anterior extent of the structure.  This
    % 'orient' flag would indicate to Analyze that this data should be placed in
    % the X-Y plane of the 3D Analyze Coordinate System, with the Z dimension
    % being the slice direction.
    
    % For the 'transverse unflipped' type, the voxels are stored with
    % Pixels in 'x' axis (varies fastest) - from patient right to left
    % Rows in   'y' axis                  - from patient posterior to anterior
    % Slices in 'z' axis                  - from patient inferior to superior
    
    fprintf('...reading axial unflipped orientation\n');
    
    avw.img = zeros(PixelDim,RowDim,SliceDim);
    
    n = 1;
    x = 1:PixelDim;
    for z = 1:SliceDim,
      for y = 1:RowDim,
        % load Y row of X values into Z slice avw.img
        avw.img(x,y,z) = tmp(n:n+(PixelDim-1));
        n = n + PixelDim;
      end
    end
    
    % no need to rearrange avw.hdr.dime.dim or avw.hdr.dime.pixdim
    
    
  case 1, % coronal unflipped
    
    % orient = 1:  The primary orientation of the data on disk is in the coronal
    % plane relative to the object scanned.  Most commonly, the fastest moving
    % index through the voxels that are part of this coronal image would span the
    % right-left extent of the structure imaged, with the next fastest moving
    % index spanning the inferior-superior extent of the structure.  This 'orient'
    % flag would indicate to Analyze that this data should be placed in the X-Z
    % plane of the 3D Analyze Coordinate System, with the Y dimension being the
    % slice direction.
    
    % For the 'coronal unflipped' type, the voxels are stored with
    % Pixels in 'x' axis (varies fastest) - from patient right to left
    % Rows in   'z' axis                  - from patient inferior to superior
    % Slices in 'y' axis                  - from patient posterior to anterior
    
    fprintf('...reading coronal unflipped orientation\n');
    
    avw.img = zeros(PixelDim,SliceDim,RowDim);
    
    n = 1;
    x = 1:PixelDim;
    for y = 1:SliceDim,
      for z = 1:RowDim,
        % load Z row of X values into Y slice avw.img
        avw.img(x,y,z) = tmp(n:n+(PixelDim-1));
        n = n + PixelDim;
      end
    end
    
    % rearrange avw.hdr.dime.dim or avw.hdr.dime.pixdim
    avw.hdr.dime.dim(2:4) = int16([PixelDim,SliceDim,RowDim]);
    avw.hdr.dime.pixdim(2:4) = single([PixelSz,SliceSz,RowSz]);
    
    
  case 2, % sagittal unflipped
    
    % orient = 2:  The primary orientation of the data on disk is in the sagittal
    % plane relative to the object scanned.  Most commonly, the fastest moving
    % index through the voxels that are part of this sagittal image would span the
    % posterior-anterior extent of the structure imaged, with the next fastest
    % moving index spanning the inferior-superior extent of the structure.  This
    % 'orient' flag would indicate to Analyze that this data should be placed in
    % the Y-Z plane of the 3D Analyze Coordinate System, with the X dimension
    % being the slice direction.
    
    % For the 'sagittal unflipped' type, the voxels are stored with
    % Pixels in 'y' axis (varies fastest) - from patient posterior to anterior
    % Rows in   'z' axis                  - from patient inferior to superior
    % Slices in 'x' axis                  - from patient right to left
    
    fprintf('...reading sagittal unflipped orientation\n');
    
    avw.img = zeros(SliceDim,PixelDim,RowDim);
    
    n = 1;
    y = 1:PixelDim;         % posterior to anterior (fastest)
    
    for x = 1:SliceDim,     % right to left (slowest)
      for z = 1:RowDim,   % inferior to superior
        
        % load Z row of Y values into X slice avw.img
        avw.img(x,y,z) = tmp(n:n+(PixelDim-1));