function [Cube,Header] = read_ctf_mri(MRIfile,MACHINEFORMAT);
%READ_CTF_MRI - Read CTF MRI format (Version 2.1)
% function [Cube,Header] = read_ctf_mri(MRIfile,MACHINEFORMAT);
%function [Cube,Header] = read_ctf_mri(MRIfile);
% MRIfile is a valid /mri file name
% Cube is a 256x256xN volume of MR scans
% Header reproduced the Header MR structure from CTF file format
% 
% function [Cube,Header] = read_ctf_mri(MRIfile,MACHINEFORMAT);
% MACHINEFORMAT is an optional string that specifies the machine format for data storage
% 'l' is for little endian
% 'b' stands for big endian (Default)
%         
% CTF MRI File Version 2.1
% The CTF MRI File format used by MRIViewer consists of a binary file with a
% 1,028 byte header. The MRI data can be in 8-bit (unsigned character) or 16-bit
% (unsigned short integer) format and consists of 256 x 256 pixel slices, stored as
% 256 contiguous sagittal slices from left to right (or right to left if head orientation
% is "left-on-right"). Each slice is stored as individual pixels starting at the
% top left corner and scanning downwards row by row. Therefore the coronal
% position is fastest changing, axial position second fastest changing and sagittal
% position slowest changing value in the file, always in the positive direction for
% each axis (see section on Head Coordinate System for axis definitions). By
% default CTF MRI files have the file extension ".mri"

%<autobegin> ---------------------- 26-May-2004 11:34:13 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Utility - General
%
% At Check-in: $Author: Mosher $  $Revision: 6 $  $Date: 5/26/04 10:02a $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 24-May-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 26-May-2004 11:34:13 -----------------------


% MRIViewer uses these cardinal directions as axes in an internal coordinate system
% where sagittal = X, coronal = Y and axial = Z forming an additional
% right-handed coordinate system which is translated and rotated with respect to
% the Head Coordinate System and has its origin at the upper left anterior corner
% of the volume.


% /---Script Authors-------------------------------------\
% |  *** Sylvain Baillet, Ph.D.                          |
% |  Cognitive Neuroscience & Brain Imaging Laboratory   |
% |  CNRS UPR640 - LENA                                  | 
% |  Hopital de la Salpetriere, Paris, France            |
% |  sylvain.baillet@chups.jussieu.fr                    |
% \------------------------------------------------------/
%  
% Script History ---------------------------------------------------------------
% SB  11 Mar-2004 Creation 
% Script History ---------------------------------------------------------------

if nargin == 1
    MACHINEFORMAT = 'b';
end


file = fopen(MRIfile,'rb',MACHINEFORMAT); 
if file < 1
    errordlg(['Error when opening : ',fullfile(pwd,MRIfile)]);
    return
end

Header.identifierString = char(fread(file,32,'char')');
Header.imageSize = fread(file,1,'short');
Header.dataSize = fread(file,1,'short');
Header.clippingRange = fread(file,1,'short');
Header.imageOrientation= fread(file,1,'short');
Header.mmPerPixel_sagittal = fread(file,1,'float32');
Header.mmPerPixel_coronal = fread(file,1,'float32');
Header.mmPerPixel_axial = fread(file,1,'float32');

% Head Model Info
Header.HeadModel_Info.Nasion_Sag = fread(file,1,'short');
Header.HeadModel_Info.Nasion_Cor = fread(file,1,'short');
Header.HeadModel_Info.Nasion_Axi = fread(file,1,'short');
Header.HeadModel_Info.LeftEar_Sag = fread(file,1,'short');
Header.HeadModel_Info.LeftEar_Cor = fread(file,1,'short');
Header.HeadModel_Info.LeftEar_Axi = fread(file,1,'short');
Header.HeadModel_Info.RightEar_Sag = fread(file,1,'short');
Header.HeadModel_Info.RightEar_Cor = fread(file,1,'short');
Header.HeadModel_Info.RightEar_Axi = fread(file,1,'short');
Header.HeadModel_Info.defaultSphereX= fread(file,1,'float32');
Header.HeadModel_Info.defaultSphereY= fread(file,1,'float32');
Header.HeadModel_Info.defaultSphereZ= fread(file,1,'float32');
Header.HeadModel_Info.defaultSphereRadius= fread(file,1,'float32');
% Head Model Info ********************************************

% Image Info
Header.Image_Info.modality = fread(file,1,'short');
Header.Image_Info.manufacturerName = fread(file,64,'char');
Header.Image_Info.instituteName = fread(file,64,'char');
Header.Image_Info.patientID = fread(file,32,'char');
Header.Image_Info.dateAndTime = fread(file,32,'char');
Header.Image_Info.scanType = fread(file,32,'char');
Header.Image_Info.contrastAgent = fread(file,32,'char');
Header.Image_Info.imagedNucleus = fread(file,32,'char');
Header.Image_Info.Frequency = fread(file,1,'float32');
Header.Image_Info.FieldStrengh = fread(file,1,'float32');
Header.Image_Info.EchoTime = fread(file,1,'float32');
Header.Image_Info.RepetitionTime = fread(file,1,'float32');
Header.Image_Info.InversionTime = fread(file,1,'float32');
Header.Image_Info.FlipAngle = fread(file,1,'float32');
Header.Image_Info.NoExcitations = fread(file,1,'short');
Header.Image_Info.NoAcquisitions = fread(file,1,'short');
Header.Image_Info.commentString = fread(file,256,'char');
Header.Image_Info.forFutureUse = fread(file,64,'char');
% Image Info ********************************************


Header.headOrigin_sagittal = fread(file,1,'float32');
Header.headOrigin_coronal = fread(file,1,'float32');
Header.headOrigin_axial  = fread(file,1,'float32');

Header.rotate_coronal = fread(file,1,'float32');
Header.headOrigin_sagittal = fread(file,1,'float32');
Header.headOrigin_axial  = fread(file,1,'float32');
Header.orthogonalFlag = fread(file,1,'short');
Header.unused = fread(file,272,'uchar');

fseek(file,1028,-1); % Skip the whole header by padding to 1028

% Header Info  ********************************************

%Read Image

tmpString = {'uchar','ushort'};

try 
    dataClass = tmpString{Header.dataSize};
catch
    errordlg('Unrecognized data format - Please check the file format');
    return
end

Cube = int16(fread(file, prod([256 256 Header.imageSize]),dataClass)); % Read the whole volume / store in a huge vector   
fclose(file);
Cube = reshape(Cube,256,256,Header.imageSize); % reshape it into a 3-D array.

