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

% function to read ESRF data format (edf & ehf)

%

% Michael Sztucki, December 2005, November 2008

%                  May 2010, Jan-Sept 2011, Feb 2012, Sept 2012

%                  July 2014

% V3.3 - with support for compressed files (zlib)

%        also uint images can be read

%

% input:     filename

% output:    image_data

%            variance_data

%            header

%

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

function [image_data,variance_data,header]=read_edf33(varargin)

if nargin == 1      

   filename=varargin{1};

   datablock=1;

   entry=-1;

   pheader='';

elseif nargin == 2

   filename=varargin{1};

   datablock=varargin{2};

   entry=-1;

   pheader='';

elseif nargin == 3

   filename=varargin{1};

   datablock=varargin{2};

   entry=varargin{3};

   pheader='';

elseif nargin == 4

   filename=varargin{1};

   datablock=varargin{2};

   entry=varargin{3};

   pheader=varargin{4};

else

   return

end

if(exist(filename)~=0)     % if the file exists

   [pathstr, name, ext] = fileparts(filename);

   variance_data=0;

   % include zipped edf images

   originalfilename=filename;

   deletetemp=0;

   deletetemp2=0;

   if strcmp(ext,'.gz') | strcmp(ext,'.GZ')

       tempfilename=[datestr(now,30) '_unzipped.edf'];

       if isunix

           unix(['gunzip -c "' filename '" > ' tempfilename]);            

       else

           unix(['"gunzip.exe" -c "' filename '" > ' tempfilename]);

       end

       filename=tempfilename;

       deletetemp=1;

   end

   fid=fopen(filename,'r');

   bo=0; % number of required header lines...

   line=' '; z=0; ehffile=0;

   while line~='{'        % find start...

       line=fgetl(fid);    

       if isempty(line) line=' '; end      % avoid warning

       z=z+1;

   end      

   line=fgetl(fid);       % read line

   header.offset_1='0';

   header.offset_2='0';

   header.center_1='0';

   header.center_2='0';

   header.detectorrotation_1='0';

   header.detectorrotation_2='0';

   header.detectorrotation_3='0';

   header.bsize_1='1';

   header.bsize_2='1';

   header.psize_1='5e-5';

   header.psize_2='5e-5';

   header.sampledistance='1';

   header.wavelength='1e-10';

   header.filename=originalfilename;

   header.edf_datablocks='0';

   header.edf_headersize='0';

   header.edf_binarysize='0';

   header.edf_binaryfileposition='0';

   header.title='';

   header.dummy='-10';    

   header.ddummy='0';

   header.projectiontype='Saxs';

   header.rasterorientation='1';

   header.axistype_1='Distance'; % can be 'Distance' (default), 'Angle', or 'Numerator'

   header.axistype_2='Distance';

   header.compression='none';

   header.edf_binaryfilename='';

   while ~(size(findstr(line,sprintf('}\n')),2)==1)

       [variable,value]=analyse_edfline(line);

       if size(variable)>0

           if isvarname(variable)

           eval(['header.' variable '=''' value ''';']);

           end

           % count required header lines

           bo=bo+strcmp(variable,'byteorder')+strcmp(variable,'datatype')+strcmp(variable,'dim_1')+strcmp(variable,'dim_2');

       end

       line=fgets(fid);

       if size(line)==1

           break

       end

   end

   edf_datablocks=0;

   if bo>=4

       if strcmp(lower(header.byteorder),'highbytefirst')

           byteorder='b';

       else

           byteorder='l';

       end

   else

       image_data=0; variance_data=0; header.title='error';fclose(fid); return;

   end

   header.edf_headersize=num2str(ftell(fid)-(size(line,2)-findstr(line,sprintf('}\n'))-1));

   edfmulti=0;

   while ~feof(fid)  && (str2num(header.edf_binarysize)>1 )

           % this might be a multi-edf file - let's search for further

           % headers...

           while ~feof(fid)

               edf_datablocks=edf_datablocks+1;

               fseek(fid,edf_datablocks*(str2num(header.edf_headersize)+str2num(header.edf_binarysize)),-1);

               line=fgets(fid);

               while line~='{'       % find start...

                   line=fgetl(fid);    

                   if line==-1

                       break

                   end

               end

               edfmulti=1;

           end

           header.edf_datablocks=num2str(edf_datablocks);

   end

   if ( edfmulti==1 ) && edf_datablocks>1

       header.title=[header.title ' (' num2str(datablock) ' of ' num2str(header.edf_datablocks) ')'];

   end

   fclose(fid);

   fid=fopen(filename,'r',byteorder); % re-open with good format

   xsize=str2num(header.dim_1);

   ysize=str2num(header.dim_2);

   rows=1:xsize;

   columns=1:ysize;

   switch lower(header.datatype)

       case {'unsignedbyte','unsigned8'}

           datatype='uint8';

           nbytes=1;

       case {'signedbyte','signed8'}

           datatype='int8';

           nbytes=1;

       case {'unsignedshort','unsigned16'}

           datatype='uint16';

           nbytes=2;

       case {'signedshort','signed16'}

           datatype='int16';

           nbytes=2;

       case {'unsignedinteger','unsignedlong','unsigned32'}

           datatype='uint32';

           nbytes=4;

       case {'signedinteger','signedlong','signed32'}

           datatype='int32';

           nbytes=4;

       case 'unsigned64'

           datatype='uint64';

           nbytes=8;

       case 'signed64'

           datatype='int64';

           nbytes=8;

       case {'float','floatvalue','real'}

           datatype='float32';

           nbytes=4;

       case {'doublevalue','double'}

           datatype='double';

           nbytes=8;

   end

   if edfmulti==1 && edf_datablocks>1

       fseek(fid,str2num(header.edf_headersize)*datablock+str2num(header.edf_binarysize)*(datablock-1),-1);               % re-position at end of header

   else

       fseek(fid,str2num(header.edf_headersize),-1);               % re-position at end of 1st header

   end

   if strcmpi(header.compression,'zcompression')

       if strcmpi(datatype,'float32')

           datatype='single';

       end

       Z=fread(fid,str2num(header.edf_binarysize));

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

       bytes=typecast(c.toByteArray,['uint8']);

       if strcmp(byteorder,'l')

           image_data=reshape((typecast(bytes,datatype)),str2num(header.dim_1),str2num(header.dim_2));

       else

           image_data=reshape(swapbytes(typecast(bytes,datatype)),str2num(header.dim_1),str2num(header.dim_2));

       end

   else

       image_data=fread(fid,[length(rows),length(columns)],sprintf('%d*%s',length(rows),datatype),nbytes*(xsize-length(rows)));

   end

   image_data=double(image_data);

   if strcmp(header.edf_binarysize,'0') || datablock~=1

       variance_data=0;

   else

       fseek(fid,str2num(header.edf_headersize)*2+str2num(header.edf_binarysize),-1);               % re-position at end of 2nd header

       if strcmpi(header.compression,'zcompression')

           if strcmpi(datatype,'float32')

               datatype='single';

           end

           Z=fread(fid); % there should be a lenghts specified

           if size(Z,1)>1

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

               bytes=typecast(c.toByteArray,'uint8');

               if strcmp(byteorder,'l')

                   variance_data=reshape((typecast(bytes,datatype)),str2num(header.dim_1),str2num(header.dim_2));

               else

                   variance_data=reshape(swapbytes(typecast(bytes,datatype)),str2num(header.dim_1),str2num(header.dim_2));

               end  

           else

               variance_data=0;

           end

       else

           variance_data=fread(fid,[length(rows),length(columns)],sprintf('%d*%s',length(rows),datatype),nbytes*(xsize-length(rows)));

       end

   end

   if max(size(image_data)~=size(variance_data))

       variance_data=0;

   end

   fclose(fid);

   if deletetemp==1

       delete(tempfilename);

   end

   if deletetemp2==1

       delete(tempfilename2);

   end

end; %end if file exists    

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

function [result1,result2]=analyse_edfline(data)

data=deblank(data);

pos=findstr('=',data);

if ~isempty(pos)

   result1=lower(strtrim(strrep(strrep(data(1:pos(1)-1),'-','_'),'~','_')));

   result2=strtrim(data(pos(1)+1:size(data,2)-1));

else

   result1='';

   result2='';

end

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