From a53aa5b3a4e1ba2d7451462f1a66559f58e12fa8 Mon Sep 17 00:00:00 2001
From: vtvivian <27768054+vtvivian@users.noreply.github.com>
Date: Thu, 26 Dec 2019 13:02:33 +0000
Subject: [PATCH] additional functionality for gridify.m

added squarify2 function, which allows step size to be defined explicitly (ebsdGrid.dx) instead of auto calculated
this sometimes leads to seam artefacts in ebsdGrid output, so all elements in ebsdgrid are assigned the nearest ebsd datapoint to remove these seams, with no interpolation of orientation data.
---
 EBSDAnalysis/@EBSD/gridify.m | 86 +++++++++++++++++++++++++++++++++++-
 1 file changed, 84 insertions(+), 2 deletions(-)

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