% read_dicom_rawdata -  Read MR raw data store in DICOM image 
%
% Marquette University,   Milwaukee, WI  USA
% Copyright 2009, 2018 - All rights reserved.
% Author: Fred J. Frigo
% Date:  May 17, 2018
%
% March 10, 2009   - Original
%                 Raw data (from a Pfile) can be added to the image data 
%                 base using the Raw File Manager "Disk to Data Base" 
%                 utility.  The new image will be created in the same exam.
%                 (e.g.  raw data from series 3 becomes series 300)
%
%                 The 'list_select_ex' utility can be used to extract
%                 this raw data image from the image data base.
%
%                 The image data base has limits for the largest image
%                 that can be stored, so very large Pfiles are not supported.
%  
% Note: the file "gems-dicom-dict.txt" (MR) can be found in this directory:
%       /export/home/sdc/app-defaults/dicom/gems-dicom-dict.txt
%

function  read_dicom_rawdata( dfile1, dict_file )

   i = sqrt(-1.0);

   if(nargin == 0)
       [fname, pname] = uigetfile('*.*', 'Select raw data image');
       dfile1 = strcat(pname, fname);
 
       % Assume DICOM dictionary file is in SAME directory for now
       dict_file = strcat(pname, 'gems-dicom-dict.txt');   
   end

   % From the DICOM filename create a NEW filename to hold raw data
   pfilename = strcat(dfile1, '.pfile');
   
   % Set dictionary to the new DICOM DICTIONARY
   dicomdict('set', dict_file); 
   dictionary = dicomdict('get');
   
   
   % Get DICOM info from image.
   info1 = dicominfo(dfile1);
   info1.UseDictionaryVR = true;
   exam = info1.StudyID;
   series = info1.SeriesNumber;
   
   image_number1 = info1.InstanceNumber;
   
   % Read Pfile header structures
   rdb_hdr_rec = info1.RdbHdrRec;
   rdb_hdr_per_pass = info1.RdbHdrPerPassTab;
   rdb_hdr_unlock_raw = info1.RdbHdrUnlockRaw;
   rdb_hdr_data_acq_tab = info1.RdbHdrDataAcqTab;
   rdb_hdr_nex_tab = info1.RdbHdrNexTab;
   rdb_hdr_nex_abort_tab = info1.RdbHdrNexAbortTab;
   rdb_hdr_tool = info1.RdbHdrTool;
   % Really need the PRESCAN_HEADER structure (14.0 and later!)
   
   %Read raw Data from the image
   raw_data = info1.RdbRawData;
   
   % Determine the Pfile Revision Number.
   f_hdr_value = typecast(info1.RdbHdrRec, 'single');
   i_hdr_value = typecast(info1.RdbHdrRec, 'uint32');
   rdbm_rev_num = f_hdr_value(1);
   if (rdbm_rev_num == 9.0)  % 11.0 product release
       pfile_header_size= 61464;
   elseif (rdbm_rev_num == 11.0)  % 12.0 product release
       pfile_header_size= 66072;
   elseif (rdbm_rev_num > 11.0) & (rdbm_rev_num < 26.0)  
       % For 14.0 to 25.0 Pfile header size can be found here
       pfile_header_size = i_hdr_value(368);   
   elseif ( rdbm_rev_num >= 26.0) & (rdbm_rev_num < 100.0)  % 26.0 to ??
       % For 26.0 Pfile header size moved to location just after rev num
       pfile_header_size = i_hdr_value(2);
   else
       err_msg = sprintf('Invalid Pfile header revision: %f', rdbm_rev_num )
       return;
   end
   
   % Create new Pfile (with PRESCAN_HEADER, EXAM, SERIES, MRIMAGE structs blank) 
   fid = fopen(pfilename,'w', 'ieee-le');
   if fid == -1
      err_msg = sprintf('Unable to open Pfile %s', pfilename)
      return;
   end
   
   % Write a blank header first since there are some missing structs
   dummy_pfile_header = zeros(1,pfile_header_size, 'uint8');                 
   status = fwrite( fid, dummy_pfile_header, 'uint8');
   
   % Create the Pfile header writing out the Pfile structures
   status = fseek(fid, 0, 'bof');
   status = fwrite(fid, rdb_hdr_rec, 'uint8');
   status = fwrite(fid, rdb_hdr_per_pass, 'uint8');
   status = fwrite(fid, rdb_hdr_unlock_raw, 'uint8');
   status = fwrite(fid, rdb_hdr_data_acq_tab, 'uint8'); 
   status = fwrite(fid, rdb_hdr_nex_tab, 'uint8'); 
   status = fwrite(fid, rdb_hdr_nex_abort_tab, 'uint8');
   status = fwrite(fid, rdb_hdr_tool, 'uint8');
   % There is a gap here at the end of the Pfile header
   % The missing data can be found in misc locations in the DICOM image header.

   % Now write the raw data.
   status = fseek(fid, 0, 'eof');
   status = fwrite(fid, raw_data, 'uint8');
   fclose(fid);
   
   msg = sprintf('New Pfile created = %s', pfilename);
   msg
  
end