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