function [ImageData, ImageDimension, ImagePrecision, VoxelDimansion] = Read_Analyze(ImagePath)
%
%Read_Analyze
%
%Reading an Analyze Image Format
%Format [ImageData, ImageSize] = UT_AnalyzeRead(ImagePath)
%
%ImagePath: The Image Full Path (It must not contain the file extension)
%ImageData: 3D Matrix Containing Image Data
%ImageDimension: Volume Dimension [X Y Z T]
%ImagePrecision: Image Data Type

if ischar (ImagePath)

    [Path, Name, Ext] = fileparts(ImagePath);
    Ext = upper (Ext);
    if (eq_str(Ext,''))||(eq_str(Ext,'.IMG'))||(eq_str(Ext,'.HDR'))
        ImagePath = strcat(Path,'/',Name,'.img');
         
    else
        ImagePath = strcat(Path,'/',Name, Ext, '.img');
    end

    if ~exist(ImagePath)
        error('File could not be found! Check the file name and path.')
    end

    AnalyzeFile = fmris_read_analyze(ImagePath);
    ImageData = AnalyzeFile.data;
    ImageDimension = AnalyzeFile.dim;
    ImagePrecision = AnalyzeFile.precision;
    VoxelDimansion(1)=AnalyzeFile.hdr.pixdim(2);
    VoxelDimansion(2)=AnalyzeFile.hdr.pixdim(3);
    VoxelDimansion(3)=AnalyzeFile.hdr.pixdim(4);

else
    error('ImagePath must be a string')
end
return;

% AnalyzeFile = fmris_read_analyze(strcat(ImagePath,'.img'));
% ImageData = AnalyzeFile;


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

% Read in 1 4D or >1 3D 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

% Get machine format on which the file was written:

machineformat=getmachineformat(file);

if ~isempty(machineformat)

    % Read Header Information
    d=Read_Analyze_Hdr(file,machineformat);
    d.file_name=file;

    % try to get file precision if it is unknown:
    if strcmp(d.precision,'Unknown')
        d.precision=getprecision(deblank(file),machineformat,d.dim,d.global);
    end

    % 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_name,'r',machineformat);
            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 all(Z>0)&all(Z<=d.dim(3))&all(T>0)&all(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(1)~=0)|(Z(1)~=0)

                if all(Z>0)&all(Z<=d.dim(3))&all(T>0)&all(T<=d.dim(4))

                    fid = fopen(d.file_name,'r',machineformat);
                    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;
else
    errordlg('Unknown Machine Format','Read Error');
end

% if there is no slice thickness, set it to 6mm:

if d.vox(3)==0
    d.vox(3)=6;
end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function machineformat=getmachineformat(file);

% Get machine format by reading the d.hdr.dim(1) attribute and
% making sure it is 1, 2, 3 or 4.

machineformat=[];
for mf='nlbdgcas'
    fid = fopen([file(1:(length(file)-3)) 'hdr'],'r',mf);
    if fid > -1
        fseek(fid,40,'bof');
        if any(fread(fid,1,'int16')==1:4)
            machineformat=mf;
            fclose(fid);
            break
        else
            fclose(fid);
        end
    else
        errordlg('Check Header File: Existence, Permissions ?','Read Error');
        break
    end
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function precision=getprecision(file,machineformat,dim,range);

% Get precision by reading a value from the middle of the .img file and
% making sure it is within the global attribute

precisions=['int8   '
    'int16  '
    'int32  '
    'int64  '
    'uint8  '
    'uint16 '
    'uint32 '
    'uint64 '
    'single '
    'float32'
    'double '
    'float64' ];
nbytes=[1 2 4 8 1 2 4 8 4 4 8 8];
middle_vol=dim(1)*dim(2)*floor(dim(3)/2)+dim(1)*round(dim(2)/2)+round(dim(1)/2);
h=dir(file);
n=ceil(h.bytes/prod(dim));

fid = fopen(file,'r',machineformat);
if fid > -1
    for i=1:size(precisions,1)
        if nbytes(i)==n
            status=fseek(fid,middle_vol*n,'bof');
            if status==0
                precision=deblank(precisions(i,:));
                x=fread(fid,10,precision);
                if all(range(1)<=x) & all(x<=range(2))
                    return
                end
            end
        end
    end
end
errordlg('Check Header File: Existence, Permissions ?','Read Error');
precision='Unknown';

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function d=Read_Analyze_Hdr(file,machineformat);

% Read Analyze Header information into the structure d
% Adapted from John Ashburners spm_hread.m

fid  = fopen([file(1:(length(file)-3)) 'hdr'],'r',machineformat);

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)));

    if ~isnan(str2double(d.hdr.scannum));    d.df(1)  =str2double(d.hdr.scannum);    end
    if ~isnan(str2double(d.hdr.patient_id)); d.df(2)  =str2double(d.hdr.patient_id); end
    if ~isnan(str2double(d.hdr.exp_date));   d.fwhm(1)=str2double(d.hdr.exp_date);   end
    if ~isnan(str2double(d.hdr.exp_time));   d.fwhm(2)=str2double(d.hdr.exp_time);   end
    if ~isnan(str2double(d.hdr.hist_un0));   d.nconj  =str2double(d.hdr.hist_un0);   end

    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

function ans = eq_str (s1,s2);
%
%		ans = eq_str (s1,s2);
%
%	returns 1 if strings are equal
%          0 if they're not
%
%		Ziad Saad July 20 97

ls1 = length (s1);
ls2 = length (s2);

if (ls1 ~= ls2),
    ans = 0;
    return;
else
    ans = 1;
    for i=1:1:ls1,
        if (s1(i) ~= s2(i)),
            ans=0;
        end;
    end;
    return;
end;


