function [ MAP, MLINE, ROIMAP ] = otl2map(X, Y, SZ, NFILL)
% [ MAP, MLINE, ROIMAP ] = otl2map(X, Y, SZ, NFILL)
%
% Try to guess which part of an image is "inside" of a (possibly complicated)
% set of contours. 
%
% NFILL (scalar) - outline points may be connected if they are NFILL+1 or fewer
%      pixels apart.
% SZ (2-vector)  - size of the image, i.e. the vector [NY,NX]
% X, Y (vectors) - x and y coordinates of all points in the outline
%
% MAP (full matrix of size SZ, and of type uint8) is a coded matrix identifying 
%   locations of the outline, and guessed locations of the inside and outside of the 
%   regions.
%              {  2 if y,x is inside the region
%   MAP(y,x) = {  1 if y,x is an outline point
%              {  0 if y,x is outside the region
%
% ROIMAP (sparse matrix of size SZ) is a matrix just like the MAP returned by otl2mline, 
%   except that each multiline object is converted to a polygon by connecting the tail 
%   point back to the head point.
%
% Algorithm:
%   1) MAP is initialized by calling otl2ml, which tries to create continuous outline 
%   contours.
%   2) Then MAP is scanned from right to left, and every time we cross
%   a contiguous series of outline points, the sign of the succeeding points
%   is toggled -- unless this contiguous series of points is missing extensions either 
%   below or above, in which case we assume that this is the "cusp" of a curve, so we 
%   shouldn't toggle.
%   3) If MAP is still +2 by the time we reach the right side, there must be an error.
%   We try to correct it by looking for the nearest scan lines above and below, and 
%   letting them vote.
%


% Initialize MAP by calling otl2ml
[MAP, MLINE] = otl2ml(X, Y, SZ, NFILL);

% Initialize the ROIMAP, then if X is zero length, return it as-is
ROIMAP = sparse(MAP);
if length(X)==0, return; end;

% Create a "troubleflags" vector to keep track of trouble rows
troubleflags = zeros(SZ(1),1);

% Create a copy of MAP in which multiline objects are replaced by closed polygons
% Call this ROIMAP.  It is made by filling pixels on the line between begin and end of each MLINE.
for i=1:length(MLINE),
  msz = size(MLINE{i});
  bpt = MLINE{i}(1,:);
  ept = MLINE{i}(msz(1),:);
  bedist = norm(bpt-ept);
  for step=1:bedist,
    pt = round((step*bpt + (bedist-step)*ept) / bedist);
    ROIMAP(pt(1), pt(2)) = 1;
  end;
end;

% Find all of the nonzero elements of ROIMAP
[YOT, XOT] = find(ROIMAP);
rows = unique(YOT);

% Since all contours are closed in ROIMAP, we know that rows 1 and SZ(1) are off-limits
rows = rows(rows~=1 & rows~=SZ(1));

% Convert MAP to non-sparse, since the 2's we add in the outlined ROI will increase the size a lot
MAP = uint8(full(MAP));

% Scan each row which has entries in ROIMAP
for row=rows(:)',
  mode = 0;
  which = find(YOT == row);
  otls = [ XOT(which); SZ(2)+1 ];

  thiso = 1;
  while thiso <= length(which),
    thisx = otls(thiso);
    nexto = thiso+1;
    nextx = thisx+1;
    
    % First, find any contiguous block of outline points on the same line
    while nexto <= length(which) & otls(nexto) <= nextx,
      nexto = nexto+1;
      nextx = nextx+1;
    end;

    % Look for extensions in ROIMAP in both up and down directions
    eup = find(ROIMAP(row-1,max(1,(thisx-1)):min(SZ(2),nextx)));
    edn = find(ROIMAP(row+1,max(1,(thisx-1)):min(SZ(2),nextx)));
    
    % Switch mode only if we have extensions in both directions
    if length(eup)>0 & length(edn)>0,
      mode = 2-mode;
    end;
    
    if mode > 0,
      % If mode is on, color this segment
      MAP(row,nextx:(otls(nexto)-1)) = 2;
      
      % If mode is on and we are at the end of the row, there is troubleflags here; mark it
      if nexto > length(which),
	troubleflags(row) = 1;
      end;
    end;

    % Finally, update thiso!
    thiso = nexto;
  end;  
end;

% Rows which were mode=2 at right-hand-side probably had errors, so we go back and fix each
trouble = find(troubleflags);
rs = 1;
while rs <= length(trouble),
  re = rs;
  while re<length(trouble) & trouble(re+1)<=trouble(re)+1,
    re = re+1;
  end;
  
  % Let rows rs-1 and re+1 "vote" to decide what to do on intermediate rows
  % We round MAP-1, so that 0+1 rounds to 0, 1+2 rounds to 2.
  % 0 on one side and 2 on the other installs an "outline point;" is that OK?
  vote = 1 + round((double(MAP(trouble(rs)-1,:)) + double(MAP(trouble(re)+1,:)) - 2)/2);
  
  for r=trouble(rs:re)',
    % Replace all elements with vote, then go back and re-fill the outline
    MAP(r,:) = uint8(vote);
    MAP(r,XOT(YOT==r)) = 1;
  end;
  
  % Finally, update rs
  rs = re+1;
end;


