function DMAP = mkdmap( xOTL, yOTL, MAP )
% DMAP = mkdmap( xOTL, yOTL, MAP )
%
% Create a DMAP of size SZ which shows how far each pixel is to the closest point in (xOTL, yOTL).
%
% xOTL - vector containing X coordinates of every point on the target outline.
% yOTL - vector containing Y coordinates of every point on the target outline.
% MAP - a MAP of type uint8, as returned by otl2map
%
% Initialize the DMAP
SZ = size(MAP);
DMAP = zeros(SZ);
% Find the rectangle which contains the OTL
rOTL = [ min(yOTL) min(xOTL) max(yOTL) max(xOTL) ];
% Find the range of possible X and Y coordinates
yRNG=rOTL(1):rOTL(3);
xRNG=rOTL(2):rOTL(4);
lx = length(xRNG);
% Find the sign of items inside the OTL rectangle
signMAP = double(MAP(yRNG, xRNG)) - 1;
% Step through the rows of the OTL rectangle, finding the minimum distance to each point from OTL
for y=yRNG,
dmap = min(repmat((y-yOTL).^2, [1 lx]) + ...
( repmat(xRNG, [length(yOTL) 1]) - repmat(xOTL, [1 lx]) ).^2 );
DMAP(y, xRNG) = signMAP(y-rOTL(1)+1,:) .* sqrt(dmap);
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now, for each cardinal direction, we expand the front rank of rOTL
% in that direction until we hit a wall.
%
% The OTL-point which is closest to any outside-point on the K'th expand must have also
% been the closest OTL-point to at least one outside-point on the (K-1)'st expand;
% we take advantage of this fact to minimize search time.
%
% Number of times we can expand in each direction before hitting a wall
numberofexpands = [ rOTL(1:2)-1 SZ-rOTL(3:4) ];
% Set the expanding direction to 1 (Y) or 2 (X); non-expanding direction is 3-dEXP.
for dEXP = 1:2,
dNEXP = 3-dEXP;
% Expand the rectangle in both - and + direction
for signEXP = [-1 1],
% Combine dEXP and signEXP information to see which direction we're going
dir = 1 + signEXP + dEXP;
pCLOSE = {yOTL xOTL};
% Now, expand as far as we can go in this particular direction
for expand=0:numberofexpands(dir),
% The expanded rectangle, clipped to the size of the image
rEXP = [ max(1, rOTL(1:2)-expand) min(SZ, rOTL(3:4)+expand) ];
% Figure out the range to cover in the non-expanding coordinate
nexpRNG = rEXP(dNEXP):rEXP(2+dNEXP);
% Create a list of the distance from each of the lRNG expansion points
% to each of the lCLOSE possible closest points
lRNG = length(nexpRNG);
lCLOSE = length(pCLOSE{1});
dmap = repmat((rEXP(dir)-pCLOSE{dEXP}).^2, [1 lRNG]) + ...
(repmat(nexpRNG, [lCLOSE 1]) - repmat(pCLOSE{dNEXP}, [1 lRNG])).^2;
% Find the shortest possible distance to each of the lRNG expansion points
if lCLOSE==1,
iCLOSEST = ones(size(dmap));
else,
[ dmap, iCLOSEST ] = min(dmap);
end
% Put dmap into the right place in DMAP, depending on whether this is Y or X
% The DMAP entries are always negative, since we know we're outside the ROI
if dEXP==1,
DMAP(rEXP(dir),nexpRNG) = - sqrt(dmap);
else,
DMAP(nexpRNG,rEXP(dir)) = - sqrt(dmap)';
end;
% Now, keep only the pCLOSE which were closest to some outside-point in the previous expand
which = unique(iCLOSEST);
pCLOSE{1} = pCLOSE{1}(which);
pCLOSE{2} = pCLOSE{2}(which);
end;
end;
end;