function S = mha_read_slice(info, num_slice)
% Function for reading the volume of a Insight Meta-Image (.mha, .mhd) file
% 
% volume = tk_read_volume(file-header)
%
% examples:
% 1: info = mha_read_header()
%    V = mha_read_volume(info);
%    imshow(squeeze(V(:,:,round(end/2))),[]);
%
% 2: V = mha_read_volume('test.mha');

if(~isstruct(info)), info=mha_read_header(info); end


switch(lower(info.DataFile))
    case 'local'
    otherwise
    % Seperate file
    info.Filename=fullfile(fileparts(info.Filename),info.DataFile);
end
        
% Open file
switch(info.ByteOrder(1))
    case ('true')
        fid=fopen(info.Filename,'rb','ieee-be');
    otherwise
        fid=fopen(info.Filename,'rb','ieee-le');
end

switch(lower(info.DataFile))
    case 'local'
        % Skip header
        fseek(fid,info.HeaderSize,'bof');
    otherwise
        fseek(fid,0,'bof');
end


datasize=prod(info.Dimensions)*info.BitDepth/8; % in bytes
slicesize=info.Dimensions(1)*info.Dimensions(2)*info.BitDepth/8;
OffsetToSelectedSlice=(num_slice-1)*slicesize;
fseek(fid, OffsetToSelectedSlice, 'cof');


switch(info.CompressedData(1))
    case 'f'
        % Read the Data
        switch(info.DataType)
            case 'char'
					S = int8(fread(fid,slicesize,'char')); 
            case 'uchar'
					S = uint8(fread(fid,slicesize,'uchar')); 
            case 'short'
					S = int16(fread(fid,slicesize,'short')); 
            case 'ushort'
					S = uint16(fread(fid,slicesize,'ushort')); 
            case 'int'
					S = int32(fread(fid,slicesize,'int')); 
            case 'uint'
					S = uint32(fread(fid,slicesize,'uint')); 
            case 'float'
					S = single(fread(fid,info.Dimensions(1)*info.Dimensions(2),'float'));   
            case 'double'
					S = double(fread(fid,slicesize,'double'));
        end
    case 't'
        switch(info.DataType)
            case 'char', DataType='int8';
            case 'uchar', DataType='uint8';
            case 'short', DataType='int16';
            case 'ushort', DataType='uint16';
            case 'int', DataType='int32';
            case 'uint', DataType='uint32';
            case 'float', DataType='single';
            case 'double', DataType='double';
        end
        Z  = fread(fid,inf,'uchar=>uint8');
        S = zlib_decompress(Z,DataType);

end

fclose(fid);
S = reshape(S,[info.Dimensions(1), info.Dimensions(2)]);

function M = zlib_decompress(Z,DataType)
import com.mathworks.mlwidgets.io.InterruptibleStreamCopier
a=java.io.ByteArrayInputStream(Z);
b=java.util.zip.InflaterInputStream(a);
isc = InterruptibleStreamCopier.getInterruptibleStreamCopier;
c = java.io.ByteArrayOutputStream;
isc.copyStream(b,c);
M=typecast(c.toByteArray,DataType);
