if h(n1)>h(n2)                                                       % Ensure grad(h)<=dhmax
         dh = (h(n1)-h(n2))/L(k);
         if dh>dhmax
            h(n1) = h(n2) + dhmax*L(k);
         end
      else
         dh = (h(n2)-h(n1))/L(k);
         if dh>dhmax
            h(n2) = h(n1) + dhmax*L(k);
         end
      end
   end
   if norm((h-h_old)./h,inf)<tol                                           % Test convergence
      break
   end
end

waitbar(0.75,wbar,'Triangulating quadtree');

t = MyDelaunayn(p);

waitbar(1,wbar);


%% SUB-FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function L = dist2poly(p,edgexy,lim)

% Find the minimum distance from the points in P to the polygon defined by
% the edges in EDGEXY. LIM is an optional argument that defines an upper
% bound on the distance for each point.

% Uses (something like?) a double sweep-line approach to reduce the number
% of edges that are required to be tested in order to determine the closest
% edge for each point. On average only size(EDGEXY)/4 comparisons need to
% be made for each point.

if nargin<3
   lim = [];
end
np = size(p,1);
ne = size(edgexy,1);
if isempty(lim)
   lim = inf*ones(np,1);
end

% Choose the direction with the biggest range as the "y-coordinate" for the
% test. This should ensure that the sorting is done along the best
% direction for long and skinny problems wrt either the x or y axes.
dxy = max(p)-min(p);
if dxy(1)>dxy(2)
    % Flip co-ords if x range is bigger
    p       = p(:,[2,1]);
    edgexy  = edgexy(:,[2,1,4,3]);
end

% Ensure edgexy(:,[1,2]) contains the lower y value
swap           = edgexy(:,4)<edgexy(:,2);
edgexy(swap,:) = edgexy(swap,[3,4,1,2]);

% Sort edges
[i,i]          = sort(edgexy(:,2));                                        % Sort edges by lower y value
edgexy_lower   = edgexy(i,:);
[i,i]          = sort(edgexy(:,4));                                        % Sort edges by upper y value
edgexy_upper   = edgexy(i,:);

% Mean edge y value
ymean = 0.5*( sum(sum(edgexy(:,[2,4]))) )/ne;

% Alloc output
L = zeros(np,1);

% Loop through points
tol = 100*eps*max(dxy);
for k = 1:np

   x = p(k,1);
   y = p(k,2);
   d = lim(k);

   if y<ymean

      % Loop through edges bottom up
      for j = 1:ne
         ymax = edgexy_lower(j,4);
         if ymax>=(y-d)
            ymin = edgexy_lower(j,2);
            if ymin<=(y+d)

               % Calculate the distance along the normal projection from [x,y]
               % to the jth edge
               x1 = edgexy_lower(j,1); x2mx1 = edgexy_lower(j,3)-x1;
               y1 = ymin;              y2my1 = ymax-y1;

               r = ((x-x1)*x2mx1+(y-y1)*y2my1)/(x2mx1^2+y2my1^2);
               if r>1                                                      % Limit to wall endpoints
                  r = 1;
               elseif r<0
                  r = 0;
               end

               dj = (x1+r*x2mx1-x)^2+(y1+r*y2my1-y)^2;
               if (dj<d^2) && (dj>tol)
                  d = sqrt(dj);
               end

            else
               break
            end
         end
      end

   else

      % Loop through edges top down
      for j = ne:-1:1
         ymin = edgexy_upper(j,2);
         if ymin<=(y+d)
            ymax = edgexy_upper(j,4);
            if ymax>=(y-d)

               % Calculate the distance along the normal projection from [x,y]
               % to the jth edge
               x1 = edgexy_upper(j,1); x2mx1 = edgexy_upper(j,3)-x1;
               y1 = ymin;              y2my1 = ymax-y1;

               r = ((x-x1)*x2mx1+(y-y1)*y2my1)/(x2mx1^2+y2my1^2);
               if r>1                                                      % Limit to wall endpoints
                  r = 1;
               elseif r<0
                  r = 0;
               end

               dj = (x1+r*x2mx1-x)^2+(y1+r*y2my1-y)^2;
               if (dj<d^2) && (dj>tol)
                  d = sqrt(dj);
               end

            else
               break
            end
         end
      end

   end

   L(k) = d;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = userhfun(x,y,fun,args,hmax,xymin,xymax)

% Evaluate user defined size function.

if ~isempty(fun)
   h = feval(fun,x,y,args{:});
   if size(h)~=size(x)
      error('Incorrect user defined size function. SIZE(H) must equal SIZE(X).');
   end
else
   h = inf*ones(size(x));
end
h = min(h,hmax);

% Limit to domain
out = (x>xymax(1))|(x<xymin(1))|(y>xymax(2))|(y<xymin(2));
h(out) = inf;

if any(h<=0)
   error('Incorrect user defined size function. H must be positive.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [hmax,edgeh,fun,args] = gethdata(hdata)

% Check the user defined size functions

if ~isempty(hdata)
   if ~isstruct(hdata)
      error('HDATA must be a structure');
   end
   if numel(hdata)~=1
      error('HDATA cannot be an array of structures');
   end
   fields = fieldnames(hdata);
   names = {'hmax','edgeh','fun','args'};
   for k = 1:length(fields)
      if ~any(strcmp(fields{k},names))
         error('Invalid field in HDATA');
      end
   end
   if isfield(hdata,'hmax')
      if (numel(hdata.hmax)~=1) || (hdata.hmax<=0)
         error('HDATA.HMAX must be a positive scalar');
      else
         hmax = hdata.hmax;
      end
   else
      hmax = inf;
   end
   if isfield(hdata,'edgeh')
      if (numel(hdata.edgeh)~=2*size(hdata.edgeh,1)) || any(hdata.edgeh(:)<0)
         error('HDATA.EDGEH must be a positive Kx2 array');
      else
         edgeh = hdata.edgeh;
      end
   else
      edgeh = [];
   end
   if isfield(hdata,'fun')
      if ~ischar(hdata.fun) && ~isa(hdata.fun,'function_handle')
         error('HDATA.FUN must be a function name or a function handle');
      else
         fun = hdata.fun;
      end
   else
      fun = '';
   end
   if isfield(hdata,'args')
      if ~iscell(hdata.args)
         error(['HDATA.ARGS must be a cell array of additional' ...
            'inputs for HDATA.FUN']);
      else
         args = hdata.args;
      end
   else
      args = {};
   end
else
   hmax = inf;
   edgeh = [];
   fun = '';
   args = {};
end