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;


