function	d=rpm_read_analyze(file,Z,T);

% Read in Analyze format Image Volumes
%
% Usage
%
% d=Read_Analyze(file);
% Reads in all image data and attributes into the structure d.
%
% d=Read_Analyze(file,Z,T);
% Reads in chosen planes and frames and attributes into the structure d.
% Z and T are vectors e.g. Z=1:31;T=1:12; or Z=[24 26];T=[1 3 5];
%
%
% (c) Roger Gunn & John Aston  


% Read Header Information
d=Read_Analyze_Hdr(file);

% Read Image Data

if ~isempty(d)&~strcmp(d.precision,'Unknown')
   
   if nargin==1|strcmp(d.precision,'uint1')
      
      % Read in Whole Analyze Volume
      fid = fopen([d.file_path d.file_name '.img'],'r','b');
      if fid > -1
         d.data=d.scale*reshape(fread(fid,prod(d.dim),d.precision),d.dim(1),d.dim(2),d.dim(3),d.dim(4));
         fclose(fid);
      else
         errordlg('Check Image File: Existence, Permissions ?','Read Error'); 
      end;
      
      if nargin==3
         if prod(Z>0)&prod(Z<=d.dim(3))&prod(T>0)&prod(T<=d.dim(4))
            d.data=d.data(:,:,Z,T);
            d.Z=Z;
            d.T=T;
         else
            errordlg('Incompatible Matrix Identifiers !','Read Error');  
         end
      end
      
   elseif nargin==3
     % Read in Chosen Planes and Frames
     if (T~=0)|(Z~=0)

     if prod(Z>0)&prod(Z<=d.dim(3))&prod(T>0)&prod(T<=d.dim(4))
         fid = fopen([d.file_path d.file_name '.img'],'r','b');
         if fid > -1
            d.data=zeros(d.dim(1),d.dim(2),length(Z),length(T));
            for t=1:length(T)
               for z=1:length(Z)
                  status=fseek(fid,d.hdr.byte*((T(t)-1)*prod(d.dim(1:3))+(Z(z)-1)*prod(d.dim(1:2))),'bof');
                  d.data(:,:,z,t)=d.scale*fread(fid,[d.dim(1) d.dim(2)],d.precision);
               end
            end
            d.Z=Z;
            d.T=T;
            fclose(fid);
         else
            errordlg('Check Image File: Existence, Permissions ?','Read Error'); 
         end;
      else
         errordlg('Incompatible Matrix Identifiers !','Read Error'); 
      end;

      end;
   else
      errordlg('Unusual Number of Arguments','Read Error');
   end;
else
   if strcmp(d.precision,'Unknown');errordlg('Unknown Data Type (Precision?)','Read Error');end
end;

return;



function d=Read_Analyze_Hdr(file);

% Read Analyze Header information into the structure d
% Adapted from John Ashburners spm_hread.m
if strcmp(computer,'PCWIN')

   separator = '\';

else

   separator			='/';
end

   
file_l    		= length(file);
if file(file_l-3)=='.';file=file(1:(file_l-4));file_l = length(file);end;
d.file_path		= file(1:max(findstr(file,separator)));
if isempty(d.file_path); 
   d.file_path	= ['.' separator];
   d.file_name	= file;
else
   d.file_name	= file((max(findstr(file,separator))+1):file_l);
end

fid   				= fopen([d.file_path d.file_name '.hdr'],'r','b');

if fid > -1
   
   % read (struct) header_key
   %---------------------------------------------------------------------------
   fseek(fid,0,'bof');
   
   d.hdr.sizeof_hdr 		= fread(fid,1,'int32');
   d.hdr.data_type  		= deblank(setstr(fread(fid,10,'char'))');
   d.hdr.db_name    		= deblank(setstr(fread(fid,18,'char'))');
   d.hdr.extents    		= fread(fid,1,'int32');
   d.hdr.session_error   	= fread(fid,1,'int16');
   d.hdr.regular    		= deblank(setstr(fread(fid,1,'char'))');
   d.hdr.hkey_un0    		= deblank(setstr(fread(fid,1,'char'))');
   
   
   
   % read (struct) image_dimension
   %---------------------------------------------------------------------------
   fseek(fid,40,'bof');
   
   d.hdr.dim    			= fread(fid,8,'int16');

   d.hdr.vox_units    		= deblank(setstr(fread(fid,4,'char'))');
   d.hdr.cal_units    		= deblank(setstr(fread(fid,8,'char'))');
   d.hdr.unused1			= fread(fid,1,'int16');
   d.hdr.datatype			= fread(fid,1,'int16');
   d.hdr.bitpix				= fread(fid,1,'int16');
   d.hdr.dim_un0			= fread(fid,1,'int16');
   d.hdr.pixdim				= fread(fid,8,'float');
   d.hdr.vox_offset			= fread(fid,1,'float');
   d.hdr.funused1			= fread(fid,1,'float');
   d.hdr.funused2			= fread(fid,1,'float');
   d.hdr.funused3			= fread(fid,1,'float');
   d.hdr.cal_max			= fread(fid,1,'float');
   d.hdr.cal_min			= fread(fid,1,'float');
   d.hdr.compressed			= fread(fid,1,'int32');
   d.hdr.verified			= fread(fid,1,'int32');
   d.hdr.glmax				= fread(fid,1,'int32');
   d.hdr.glmin				= fread(fid,1,'int32');
   
   % read (struct) data_history
   %---------------------------------------------------------------------------
   fseek(fid,148,'bof');
   
   d.hdr.descrip			= deblank(setstr(fread(fid,80,'char'))');
   d.hdr.aux_file			= deblank(setstr(fread(fid,24,'char'))');
   d.hdr.orient				= fread(fid,1,'char');
   d.hdr.origin				= fread(fid,5,'uint16');
   d.hdr.generated			= deblank(setstr(fread(fid,10,'char'))');
   d.hdr.scannum			= deblank(setstr(fread(fid,10,'char'))');
   d.hdr.patient_id			= deblank(setstr(fread(fid,10,'char'))');
   d.hdr.exp_date			= deblank(setstr(fread(fid,10,'char'))');
   d.hdr.exp_time			= deblank(setstr(fread(fid,10,'char'))');
   d.hdr.hist_un0			= deblank(setstr(fread(fid,3,'char'))');
   d.hdr.views				= fread(fid,1,'int32');
   d.hdr.vols_added			= fread(fid,1,'int32');
   d.hdr.start_field		= fread(fid,1,'int32');
   d.hdr.field_skip			= fread(fid,1,'int32');
   d.hdr.omax				= fread(fid,1,'int32');
   d.hdr.omin				= fread(fid,1,'int32');
   d.hdr.smax				= fread(fid,1,'int32');
   d.hdr.smin				= fread(fid,1,'int32');
   
   fclose(fid);
   
   % Put important information in main structure
   %---------------------------------------------------------------------------
   
   d.dim    	  				= d.hdr.dim(2:5)';
   vox 						= d.hdr.pixdim(2:5)';
   if 	vox(4)==0 vox(4)=[];end
   d.vox       				= vox;
   d.vox_units       		= d.hdr.vox_units;
   d.vox_offset	    		= d.hdr.vox_offset;
   scale     				= d.hdr.funused1;
   d.scale     			  	= ~scale + scale;
   d.global					= [d.hdr.glmin d.hdr.glmax];
   d.calib					= [d.hdr.cal_min d.hdr.cal_max];
   d.calib_units			= d.hdr.cal_units;
   d.origin    				= d.hdr.origin(1:3)';
   d.descrip   				= d.hdr.descrip(1:max(find(d.hdr.descrip)));
   
   switch d.hdr.datatype
   case 1
      d.precision 	= 'uint1';
      d.hdr.byte 	= 0;
   case 2
      d.precision 	= 'uint8';
      d.hdr.byte 	= 1;
   case 4
      d.precision 	= 'int16';
      d.hdr.byte 	= 2;
   case 8
      d.precision 	= 'int32';
      d.hdr.byte 	= 4;
   case 16
      d.precision 	= 'float';
      d.hdr.byte 	= 4;
   case 64
      d.precision 	= 'double';
      d.hdr.byte 	= 8;
   otherwise
      d.precision 	= 'Unknown';
      d.hdr.byte 	= 0;
   end
   
else
   d=[];
   errordlg('Check Header File: Existence, Permissions ?','Read Error'); 
end


return;