Fixed missing triangle intersections when using pointLocation

due to floating point errors.
The intersections are evaluated after shifting the reference and query values by the same offset in different directions which reduces the number of errors.

gridCoords.m

@@ -2,10 +2,34 @@ function [Xq,Yq,X,Y,varargout] = gridCoords(X,Y,Xq,Yq,varargin)
 %GRIDCOORDS Summary of this function goes here
 %   Detailed explanation goes here
 
+I = (1:numel(X))';
+offset = 1e-2;
+
+p = inputParser();
+p.addOptional('I'     , I     , @isnumeric);
+p.addOptional('offset', offset, @isnumeric);
+
+p.parse(varargin{:});
+I      = p.Results.I;
+offset = p.Results.offset;
+
 sz = size(Xq);
 
-[t,I] = utils.triangulation(X,Y,varargin{:});
-[ind,b] = t.pointLocation(Xq(:),Yq(:));
+ind = nan(numel(Xq),1);
+b   = nan(numel(Xq),3);
+nanmask = true(numel(Xq),1);
+
+[dX,dY] = ndgrid([-offset,0,offset]);
+sigma = [5,4,6,2,8,1,3,7,9];
+
+for r = 1:numel(sigma)
+    dx = dX(sigma(r));
+    dy = dY(sigma(r));
+    
+    t = utils.triangulation(X+dx,Y+dy,I);
+    [ind(nanmask),b(nanmask,:)] = t.pointLocation(Xq(nanmask)+dx,Yq(nanmask)+dy);
+    nanmask = isnan(ind);
+end
 
 ind = reshape(ind,sz);
 u = reshape(b(:,2),sz);
@@ -13,7 +37,6 @@ v = reshape(b(:,3),sz);
 
 nanind = isnan(ind);
 ind(nanind) = 1;
-ind = I(ind);
 ut = ~mod(ind,2);
 ind(nanind) = nan;
 

triangulation.m

-function varargout = triangulation(X,Y,varargin)
+function t = triangulation(X,Y,varargin)
 %TRIANGULATION Summary of this function goes here
 %   Detailed explanation goes here
 
+I = (1:numel(X))';
 p = inputParser();
-p.addOptional('Z'    , []    , @isnumeric);
-p.addOptional('order', 'none', @ischar);
-p.addOptional('fun'  , @(x) min(x,[],2));
+p.addOptional('I', I, @isnumeric);
 
 p.parse(varargin{:});
-Z = p.Results.Z;
-order = p.Results.order;
-fun = p.Results.fun;
+I = p.Results.I;
 
 sz = size(X);
 
@@ -20,26 +17,8 @@ TY = reshape( (Y_(:)+[0,0,1,1,1,0])' ,3,[])';
 
 T = sub2ind(sz,TX,TY);
 
-switch order
-    case 'none'
-        I = (1:size(T,1))';
-    otherwise
-        Z_ = Z(T);
-        Z_ = fun(Z_);
-        [~,I] = sort(Z_,order);
-        T = T(I,:);
-end
-
-if nargout>0
-    varargout{1} = triangulation(T,X(:),Y(:));
-end
-
-if nargout>1
-    varargout{2} = I;
-end
-
-if nargout>2
-    varargout{3} = triangulation(T,X(:),Y(:),Z(:));
-end
+Ii(I) = 1:numel(I);
+
+t = triangulation(Ii(T),X(I(:)),Y(I(:)));
 
 end
\ No newline at end of file