diff --git a/EBSDAnalysis/@EBSD/gridify.m b/EBSDAnalysis/@EBSD/gridify.m
index 9ea762982..9b36d4e5d 100644
--- a/EBSDAnalysis/@EBSD/gridify.m
+++ b/EBSDAnalysis/@EBSD/gridify.m
@@ -29,7 +29,9 @@
 %   plot(ebsdMg, ebsdMg.orientations)
 %
 
-if size(ebsd.unitCell,1) == 6
+if nargin==2
+    [ebsdGrid,newId] = squarify2(ebsd,varargin{:});%Vivian Dec 2019, extra function defining grid step size
+elseif size(ebsd.unitCell,1) == 6
   [ebsdGrid,newId] = hexify(ebsd,varargin{:});
 else
   [ebsdGrid,newId] = squarify(ebsd,varargin{:});
@@ -79,6 +81,86 @@
 
 end
 
+
+function [ebsdGrid,newId] = squarify2(ebsd,varargin)
+%add extra functionality to set step size and remove seam artefacts (Dec 2019)
+
+% generate regular grid
+prop = ebsd.prop;
+ext = ebsd.extend; % ebsd spatial boudnings - [xmin xmax ymin ymax]
+%assign the step size explicitly, intead of calculating from ebsd.unitCell
+dx = varargin{1};
+dy = varargin{1};%square grid
+[prop.x,prop.y] = meshgrid(linspace(ext(1),ext(2),1+round((ext(2)-ext(1))/dx)),...
+    linspace(ext(3),ext(4),1+round((ext(4)-ext(3))/dy))); % ygrid runs first
+sGrid = size(prop.x); %new grid
+
+
+% VT 2019 - this assigns every old ebsd point a location in sGrid
+% VT 2019 - which leads to unfilled seams in the final map
+
+% detect position within grid
+newId = sub2ind(sGrid, 1 + round((ebsd.prop.y - ext(3))/dy), ...
+  1 + round((ebsd.prop.x - ext(1))/dx));
+
+% set phaseId to notIndexed at all empty grid points
+phaseId = nan(sGrid);
+phaseId(newId) = ebsd.phaseId;
+
+% update rotations
+a = nan(sGrid); b = a; c = a; d = a;
+a(newId) = ebsd.rotations.a;
+b(newId) = ebsd.rotations.b;
+c(newId) = ebsd.rotations.c;
+d(newId) = ebsd.rotations.d;
+
+% update all other properties
+for fn = fieldnames(ebsd.prop).'
+    if any(strcmp(char(fn),{'x','y','z'})), continue;end
+    if isnumeric(prop.(char(fn)))
+        prop.(char(fn)) = nan(sGrid);
+    else
+        prop.(char(fn)) = prop.(char(fn)).nan(sGrid);
+    end
+    prop.(char(fn))(newId) = ebsd.prop.(char(fn));
+end
+
+%to avoid the seam artefact, assign every sGrid point its closest
+%ebsd point
+%but it's too computationally expensive to do full map this way!
+%try assigning only to previously unassigned points instead
+asg = nan(sGrid);
+asg(newId) = ebsd.id;
+[unasgR,unasgC] = find(isnan(asg)); %these are the missing points not assigned by newId
+yOld = 1 + (ebsd.prop.y - ext(3))/dy;
+xOld = 1 + (ebsd.prop.x - ext(1))/dx;
+ix = zeros(size(unasgR));
+
+for n=1:length(unasgR)
+        %to get the correct element, minimise the product of the row/colum
+        %differences
+        % min(abs(yOld - rNew))would give the nearest row
+        % min(abs(xOld - cNew)) would give the nearest column
+        % minimise sum of squares to both?
+        [~,ix(n)] = min(((yOld - unasgR(n)).^2)+((xOld - unasgC(n)).^2));
+        
+        phaseId(unasgR(n),unasgC(n)) = ebsd.phaseId(ix(n));
+        a(unasgR(n),unasgC(n)) = ebsd.rotations.a(ix(n));
+        b(unasgR(n),unasgC(n)) = ebsd.rotations.b(ix(n));
+        c(unasgR(n),unasgC(n)) = ebsd.rotations.c(ix(n));
+        d(unasgR(n),unasgC(n)) = ebsd.rotations.d(ix(n));
+        
+end
+
+ebsdGrid = EBSDsquare(rotation(quaternion(a,b,c,d)),phaseId(:),...
+    ebsd.phaseMap,ebsd.CSList,[dx,dy],'options',prop);
+
+% newId output isn't meaningful for squarify2(), but keeping it
+% as it is still consistent with squarify()
+
+end
+
+
 function [ebsdGrid,newId] = hexify(ebsd,varargin)
 
 prop = ebsd.prop;
@@ -168,4 +250,4 @@
 ebsdGrid = EBSDhex(rot,phaseId(:),...
   ebsd.phaseMap,ebsd.CSList,dHex,isRowAlignment,offset,'options',prop);
 
-end
\ No newline at end of file
+end