diff --git a/gridCoords.m b/gridCoords.m
index 69ec322eef31b32415320e3b9c78ee785f3ebb3a..f1dac7856e97bbf5c62226e08d4d691821cc655e 100644
--- a/gridCoords.m
+++ b/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;
 
diff --git a/triangulation.m b/triangulation.m
index ab50c8025887bd53e7540d2be24e87046d634886..7ce095feb8d74aa99a6188b4aa3ca3260f92bbcc 100644
--- a/triangulation.m
+++ b/triangulation.m
@@ -1,16 +1,13 @@
-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