function OK = fread_ABI_coll(gel, XL, VERBOSE)
%
% OK = fread_ABI_coll(gel, XL, VERBOSE)
%
%
% Written by See-Kiong Ng, 1997, 1998
% Copyright (c) 1997-98 by See-Kiong Ng, CMU.

% Note: To put back automatic start scan: search for commented
%       start_scan and SIZE_DATA -- uncomment them; and
%       uncomment the line "start_scan = 1;".

if nargin<2 XL = 0; VERBOSE = 0; end

ALLOW_CORRUPT = XL; 

global SYS_ARCHITECTURE 

SMOOTH = 0;



[NUM_ABI_FILTERS,DATA_PRECISION,ABI_ARCHITECTURE,MAC_HEADER_SIZE]=ABI_info;
 MAC_ARCHITECTURE = ABI_ARCHITECTURE;

if XL ABI_SCAN_WIDTH  = (1+1)*194; end    % SKNG 7/24/98
ABI_HEADER_SIZE = ABI_SCAN_WIDTH*2+17+2;
if XL==96 ABI_HEADER_SIZE = ABI_HEADER_SIZE + 184; end  % SKNG 7/24/98
ABI_RECORD_SIZE = ABI_HEADER_SIZE;


MAX_MEMORY = ABI_SCAN_WIDTH*1000;

FRONT_OFFSET = 0;

OK = 0;

colfile = path_inputfile('gel_file', gel);  % Assumes we have already 
                                            % tested is_ABI_col(gel)

[gel_machine, dye_planes, size_plane_num] = gel_overview(gel);

          
[fid, msg] = fopen(colfile, 'r', ABI_ARCHITECTURE);
if fid<0, 
   fprintf('\nError: %s - %s.\n', colfile, msg); return;
end





% STEP 1: Reads the 1st header block
% ------

header_num = 1;

% Attempts to read the first header block
%
[header, count] = fread(fid, 12, 'char');
if count<12, 
  fprintf('\nError: %s corrupted. (-1)\n', colfile); return;
end



if (header(1)~='V')|(header(12)~='F')
   % Possible Macbinary format: attempt to repair 
   fseek(fid, MAC_HEADER_SIZE, 'bof');
   [header, count] = fread(fid, 12, 'char');
   if count<12,   
     fprintf('\nError: %s corrupted. (-2)\n', colfile); return;
   end
   FRONT_OFFSET = MAC_HEADER_SIZE;
end

offset = FRONT_OFFSET;

if (header(1)~='V')|(header(12)~='F')
   fprintf(['\nError: ' colfile ' is not a valid ABI 373 collection file -- ', ...
          'corrupted header 1.\n']);
   return;
end

% Move to byte 17
%

fseek(fid, offset+17, 'bof');
[ABI_version, count] = fscanf(fid, '%s', 1);
if count<1, 
  fprintf('\nError: %s corrupted. (-3)\n', colfile); return;
end


% STEP 2: Reads the 2nd header block
% ------
%
header_num = 2;
offset = FRONT_OFFSET+(header_num-1)*ABI_HEADER_SIZE;
fseek(fid, offset, 'bof');
[header, count] = fread(fid, 12, 'char');
if count<12, 
  fprintf('\nError: %s corrupted. (-4)\n', colfile); return;
end

if (header(1)~='T')|(header(12)~='F')
   fprintf(['\nError: ' colfile ' is not a valid ABI 373 collection file -- ', ...
          'corrupted header 2.\n']);
   return;
end

% 2nd header block has no interesting datum.


% STEP 3: Reads the 3rd header block
% ------
%
header_num = 3;
offset = FRONT_OFFSET+(header_num-1)*ABI_HEADER_SIZE;
fseek(fid, offset, 'bof');
[header, count] = fread(fid, 12, 'char');
if count<12, 
  fprintf('\nError: %s corrupted. (-5)\n', colfile); return;
end
if (header(1)~='M')|(header(12)~='F')
   fprintf(['\nError: ' colfile ' is not a valid ABI 373 collection file -- ', ...
          'corrupted header 3.\n']); 
   return;
end

% 3rd header *may* contain the dye matrix file name -- ignore



% STEP 4: Reads the 4th header block
% ------
%
header_num = 4;
offset = FRONT_OFFSET+(header_num-1)*ABI_HEADER_SIZE;
fseek(fid, offset, 'bof');
[header, count] = fread(fid, 12, 'char');
if count<12, 
  fprintf('\nError: %s corrupted. (-6)\n', colfile); return;
end
if (header(1)~='C')|(header(12)~='F')
   fprintf(['\nError: ' colfile ' is not a valid ABI 373 collection file -- ', ...
          'corrupted header 4.\n']); 
   return;
end

% 4th header should point to NULL -- ignore


% STEP 5: Reads the 5th header block
% ------
%
header_num = 5;
offset = FRONT_OFFSET+(header_num-1)*ABI_HEADER_SIZE;
fseek(fid, offset, 'bof');
[header, count] = fread(fid, 12, 'char');
if count<12, 
  fprintf('\nError: %s corrupted. (-7)\n', colfile); return;
end

if (header(1)~='N')|(header(12)~='F')
   fprintf(['\nError: ' colfile ' is not a valid ABI 373 collection file -- ', ...
          'corrupted header 5.\n']); 
   return;
end

% 5th header should tell number of blocks
fseek(fid, offset+17, 'bof');
[num_ABI_blocks, count] = fread(fid, 1, 'long');
if count~=1 
  fprintf('\nError: %s corrupted. (-8)\n', colfile); return;
end

if num_ABI_blocks==2 
   NUM_HEADER_BLOCKS = 9;  % For ABI 373 A
elseif num_ABI_blocks==3
   NUM_HEADER_BLOCKS = 11; % For ABI 377
else
   fprintf(['\nError: ' colfile ' is not a valid ABI collection file -- ', ...
          'corrupted header 5 (data_block_count = ' ...
          int2str(num_ABI_blocks) ').\n']); 
  % SKNG - 3/1/98 for Genethon
  fprintf('Warning - Attempt to continue: gel file may be corrupted.\n');
  fprintf('Warning - Assuming %s is an ABI/373 gel.\n', get_filename(colfile));
  NUM_HEADER_BLOCKS = 9;

   %return;
end

% STEP 6: Reads the 6th header block
% ------
%
header_num = 6;
offset = FRONT_OFFSET+(header_num-1)*ABI_HEADER_SIZE;
fseek(fid, offset, 'bof');
[header, count] = fread(fid, 12, 'char');
if count<12, 
  fprintf('\nError: %s corrupted. (-9)\n', colfile); return;
end
if (header(1)~='S')|(header(12)~='F')
   fprintf(['\nError: ' colfile ' is not a valid ABI 373 collection file -- ', ...
          'corrupted header 6.\n']); return;
end

% 6th header tells tells the offset to primary data block
%
fseek(fid, offset+17, 'bof');
[data_offset, count] = fread(fid, 1, 'long');
if count~=1 
  fprintf('\nError: %s corrupted. (-10)\n', colfile); return;
end

if data_offset ~= (NUM_HEADER_BLOCKS*ABI_HEADER_SIZE)
   fprintf(['\nError: ' colfile ' is not a valid ABI 373 collection file -- ', ...
          'corrupted header 6 (read data offset = ' ...
          int2str(data_offset) '; expected = ' ...
          int2str((NUM_HEADER_BLOCKS*ABI_HEADER_SIZE)) ').\n']);
  % SKNG - 3/1/98 for Genethon
  fprintf('Warning - Attempt to continue: gel file may be corrupted.\n');
%  return;
end


% STEP 7: Reads the 7th header block
% ------
%
header_num = 7;
offset = FRONT_OFFSET+(header_num-1)*ABI_HEADER_SIZE;
fseek(fid, offset, 'bof');
[header, count] = fread(fid, 12, 'char');
if count<12, 
  fprintf('\nError: %s corrupted. (-11)\n', colfile); return;
end
if (header(1)~='L')|(header(12)~='F')
   fprintf(['\nError: ' colfile ' is not a valid ABI 373 collection file -- ', ...
          'corrupted header 7.\n']); return;
end

% 7th header tells tells the offset to primary data block
%
fseek(fid, offset+17, 'bof');
[data_block_len, count] = fread(fid, 1, 'long');
if count~=1 
  fprintf('\nError: %s corrupted. (-12)\n', colfile); return;
end



% STEP 8: Reads the 8th header block
% ------
%
header_num = 8;
offset = FRONT_OFFSET+(header_num-1)*ABI_HEADER_SIZE;
fseek(fid, offset, 'bof');
[header, count] = fread(fid, 12, 'char');
if count<12, 
  fprintf('\nError: %s corrupted. (-13)\n', colfile); return;
end

if (header(1)~='S')|(header(12)~='F')
   fprintf(['\nError: ' colfile ' is not a valid ABI 373 collection file -- ', ...
          'corrupted header 8.\n']); return;
end

% 8th header tells tells the offset to secondary data block
%
fseek(fid, offset+17, 'bof');
[data_offset2, count] = fread(fid, 1, 'long');
if count~=1 
  fprintf('\nError: %s corrupted. (-14)\n', colfile); return;
end
if data_offset2 ~= (data_offset+data_block_len)
   fprintf(['\nError: ' colfile ' is not a valid ABI 373 collection file -- ', ...
          'corrupted header 8 (read data offset = ' ...
          int2str(data_offset) '; expected = ' ...
          int2str((data_offset+data_block_len)) ').\n']); return;
end


% STEP 9: Reads the 9th header block
% ------
%
header_num = 9;
offset = FRONT_OFFSET+(header_num-1)*ABI_HEADER_SIZE;
fseek(fid, offset, 'bof');
[header, count] = fread(fid, 12, 'char');
if count<12, 
  fprintf('\nError: %s corrupted. (-15)\n', colfile); return;
end

if (header(1)~='L')|(header(12)~='F')
   fprintf(['\nError: ' colfile ' is not a valid ABI 373 collection file -- ', ...
          'corrupted header 9.\n']); return;
end

% 9th header tells tells the offset to secondary data block
%
fseek(fid, offset+17, 'bof');
[data_block_len2, count] = fread(fid, 1, 'long');
if count~=1 
  fprintf('\nError: %s corrupted. (-16)\n', colfile); return;
end

gel_width = ABI_SCAN_WIDTH;
num_planes = NUM_ABI_FILTERS;
gel_len = data_block_len/(num_planes*ABI_RECORD_SIZE);

if VERBOSE 
   fprintf('\nImage file %s is an ABI', get_filename(colfile));
   if XL if XL==96 fprintf('-96'); end; fprintf('/XL'); 
   elseif NUM_HEADER_BLOCKS==9, fprintf('/373');
   elseif NUM_HEADER_BLOCKS==11, fprintf('/377'); 
   else fprintf('[Unknown]'); end

   fprintf(' collection file (version %s).\n', ABI_version);
   fprintf('%s is a %i x %i gel with %i dye planes.\n\n', ...
            gel, gel_len, gel_width, num_planes);
end


mtxfile = path_inputfile('matrix_file', gel);
if isempty(mtxfile) 
   fprintf('\nError: Please specify the ABI matrix file and restart.\n');
   return;
end

if ~fexists(mtxfile) 
   fprintf('\nError: Matrix file %s does not exist.\n', mtxfile);
   return;
end

MATRIX = fread_ABI_mtx(mtxfile);
MATRIX = MATRIX/1000;

%MATRIX = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];



% STEP 10: Read each scan data block
% -------

if VERBOSE
  fprintf('Reading gel data from %s ....', get_filename(colfile));
end

data_offset = data_offset+FRONT_OFFSET;
fseek(fid, data_offset, 'bof');


fids = zeros(1, NUM_ABI_FILTERS);
start_scan = 1;  

for p=1:NUM_ABI_FILTERS
  tmpfile = [path_displayfile('plane', gel, p) '.TMP'];
  [fids(p), msg] = fopen(tmpfile, 'w', SYS_ARCHITECTURE);
  if fids(p)<0 
    fprintf('\nError: %s - %s.\n', tmpfile, msg); return;
  else
   %
   % First line contains three things:
   %    DATA_PRECISION GEL_LENGTH  GEL_WIDTH [FIRST_SCAN]
   %
%   fprintf(fids(p), '%s\t%i\t%i\t%i\n', DATA_PRECISION, ...
%           gel_len-start_scan+1, gel_width, start_scan);
%
  end
end

scan_data = zeros(NUM_ABI_FILTERS, ABI_SCAN_WIDTH);
deltas    = zeros(NUM_ABI_FILTERS, ABI_SCAN_WIDTH);
lasts     = zeros(NUM_ABI_FILTERS, ABI_SCAN_WIDTH);
sods      = zeros(NUM_ABI_FILTERS, ABI_SCAN_WIDTH);
offs      = 65000*ones(NUM_ABI_FILTERS, ABI_SCAN_WIDTH);

old_scan_data = scan_data;

% SIZE_DATA = zeros(gel_len, gel_width);

scan_line_num = 0; plane = 0; trav_count = zeros(1, NUM_ABI_FILTERS);
scan_count = 0; total_num_scans = data_block_len/ABI_RECORD_SIZE;

num_scans_at_once = ceil(MAX_MEMORY/ABI_RECORD_SIZE);

%all_record = [];

for k=1:num_scans_at_once:total_num_scans

  num_scans_now = min(num_scans_at_once, total_num_scans-k+1);
%  clear all_record
  [all_record, count] = fread(fid, ABI_RECORD_SIZE*num_scans_now, ...
                              'uchar'); 
  all_record = all_record';
  if count~=ABI_RECORD_SIZE*num_scans_now, 
     fprintf(['\nError: ' colfile ': corrupted at scan record ' ...
             int2str(scan_line_num) ...
             '(count = ' int2str(count) '; ABI_RECORD_SIZE = ' ...
             int2str(ABI_RECORD_SIZE) ').\n']); return;
  end
  if VERBOSE fprintf('.'); end

%  for scan_line_num = (k-1)+[1:num_scans_now]


  while scan_line_num<k+num_scans_now-1
     scan_line_num = scan_line_num+1;
     record = ...
      all_record((scan_line_num-k)*ABI_RECORD_SIZE+[1:ABI_RECORD_SIZE]);

     % Check some of the header info
     CORRUPTED = 0;
     if sscanf(setstr(record(2:4)), '%d',1) ~= (ABI_RECORD_SIZE-4-2)
        fprintf(['\nError: ' colfile ': primary data segment '...
                 'corrupted on scan %d.\n'], scan_line_num); 
        add_WARN(gel_or_test_indx(gel), scan_line_num, 1);
        if ~ALLOW_CORRUPT return; 
        else CORRUPTED = 1; 
          fprintf(['       - Results may be incorrect ', ...
                   'due to corrupted data.\n']);
        end
     end

     traverse_counter = sscanf(setstr(record(13:17)), '%d', 1);
   
     % Filter number ([0, 1, 2, 3])
     %

     new_plane = sscanf(setstr(record(12)), '%d', 1)+1;
     if plane==NUM_ABI_FILTERS expected_plane = 1;
     else expected_plane = plane+1; end
          
     READ_IT = 1;


     if ~plane, plane = new_plane;
     elseif new_plane~=expected_plane
       % Filter out of order


       % Mark the following filter as read.
       %
       fprintf('\nError: Filter data out of order : %d -> %d (expected %d)\n', ...
                plane, new_plane, expected_plane);
       fprintf('       - Attempt to repair by assuming the other filters read.\n');
       fprintf('       - Results may be inaccurate.\n');

       if new_plane>expected_plane
          % Skips the intermediate planes
          %
          plane = new_plane;
          trav_count(expected_plane:new_plane-1) = ...
               -999*ones(size(trav_count(expected_plane:new_plane-1)));
       else
          % Reset to next round
          %
          trav_count(expected_plane:NUM_ABI_FILTERS) = ...
               -999*ones(size(trav_count(expected_plane:NUM_ABI_FILTERS)));
          READ_IT = 0; 
          scan_line_num = scan_line_num-1; 
          plane = NUM_ABI_FILTERS;
       end

     else plane = new_plane; end

     if READ_IT
       scan_line = reshape(record([1:2*ABI_SCAN_WIDTH]+17), ...
                           2, ABI_SCAN_WIDTH);
       if rem(plane,2)==0, scan_line = fliplr(scan_line); end
       scan_line = scan_line(1,:)+scan_line(2,:)*256;
       scan_line = scan_line(1,:);
%       scan_line = scan_line - min(scan_line);

       scan_data(plane, :) = scan_line;
       trav_count(plane) = trav_count(plane)+1;
     end


     if all(trav_count~=0)
       bad_indxes = find(trav_count~=1);
       if ~isempty(bad_indxes)
           fprintf('\nWarning: Gel scan line out of order < '); 
           fprintf('%d ',trav_count);
           fprintf('>.\n');
          % Attemp to repair
           scan_data(bad_indxes,:) = old_scan_data(bad_indxes,:);
       end
       old_scan_data = scan_data;
       trav_count = zeros(1,NUM_ABI_FILTERS);
       scan_count = scan_count+1;

       % Do some smoothing (if requested)
       %
       if SMOOTH
        deltas = scan_data - lasts;
        lasts  = scan_data;
        sods   = sods + deltas;
        offs   = min(offs, sods);
        scan_data = sods - offs;
       end

       % Apply dye separation matrix
       %
       scan_data = MATRIX*scan_data;
       scan_data = scan_data.*(scan_data>0);


       %
       % And then save into planefile
       %
       for p = 1:NUM_ABI_FILTERS
         count = fwrite(fids(p), scan_data(p,:), DATA_PRECISION);
         if count~=ABI_SCAN_WIDTH
           fprintf('\nError: Failure in writing temp plane %d.\n',p);
           return;
         end
       end

       scan_data = zeros(NUM_ABI_FILTERS, ABI_SCAN_WIDTH);
     end
  end
end  % for k


for p = 1:NUM_ABI_FILTERS, fclose(fids(p)); end


if (scan_count*NUM_ABI_FILTERS ~= total_num_scans)
   fprintf('\nWarning (possible error): Scan counts = %d, but expected total %d scans.\n',...
           scan_count*NUM_ABI_FILTERS, total_num_scans);
   %return;
end

if ftell(fid) ~= data_offset+data_block_len
  fprintf(['\nError: ftell in fread_ABI_COLL - ' ...
          int2str(ftell(fid)),' ', ...
          int2str(data_offset+data_block_len) '.\n']);
       
  return;
end
fclose(fid);


%if VERBOSE fprintf('Done.\n'); end


% Now, try to determine the start scan line from the size_plane
%
%start_scan = find_ABI_start_scan(SIZE_DATA);
start_scan = 1;  

if VERBOSE fprintf('.'); end

% Finally, write the cropped image into planefile
%

for p=1:NUM_ABI_FILTERS
   planefile = path_displayfile('plane', gel, p);
   tmpfile = [planefile '.TMP'];
   [fid, msg] = fopen(tmpfile, 'r', SYS_ARCHITECTURE);
   if fid<0 
       fprintf('\nError: %s : %s.\n', tmpfile, msg);
       return;
   end

   if VERBOSE fprintf(':'); end

   NUM_ROWS_AT_ONCE = ceil(MAX_MEMORY/gel_width);
   plane = zeros(gel_len, gel_width);
   max_r = ceil(gel_len/NUM_ROWS_AT_ONCE); total_count = 0;

   for r=1:max_r
     num_rows = min(gel_len - (r-1)*NUM_ROWS_AT_ONCE, NUM_ROWS_AT_ONCE);
     [data, count] = fread(fid, [gel_width num_rows], DATA_PRECISION);
            
     if count~=num_rows*gel_width, 
        fprintf(['\nError: Corrupted %s ', ...
                 '(num_rows = %d, gel_width = %d, count = %d, ',...
                 'r = %d, max_r = %d).\n'], ...
                 get_filename(tmpfile), num_rows,gel_width,count,...
                 r, max_r); 
        return;
     else
        total_count = total_count+count;
        plane([1:num_rows]+(r-1)*NUM_ROWS_AT_ONCE,:) = data';
     end
   end

   if total_count~=gel_len*gel_width, 
     fprintf(['PROGRAM Error 986776: Corrupted ' tmpfile '.\n']); 
     return;
   end
   fclose(fid); fremove(tmpfile);

   [fid, msg] = fopen(planefile, 'w', SYS_ARCHITECTURE);
   if fid<0 
     fprintf(['\nError: ' planefile '- ' msg '.\n']);
     return;
   end

   %
   % First line contains three things:
   %    DATA_PRECISION GEL_LENGTH  GEL_WIDTH [FIRST_SCAN]
   %
   fprintf(fid, '%s\t%i\t%i\t%i\n', DATA_PRECISION, ...
           gel_len-start_scan+1, gel_width, start_scan);
   count = fwrite(fid, plane(start_scan:gel_len,:), DATA_PRECISION);
   if count ~=(gel_len-start_scan+1)*gel_width
      fprintf(['\nError: Failure in writing ' planefile '.\n']);
      return;
   end
   if VERBOSE fprintf('.'); end

   fclose(fid);          
end



if VERBOSE 
   fprintf('Done.'); 
   if start_scan>1 fprintf(' (Start scan = %d)', start_scan); end
   fprintf('\n');
end





OK = 1;
