contour based gB reconstruction and smoothing

EBSDAnalysis/@grain2d/boundaryContours.m

@@ -0,0 +1,73 @@
+function [grainsContour] = boundaryContours(grains)
+% Redraw grain boundary segments using a contouring algorithm instead of
+% following pixel edges
+%
+% This is similar to the marching squares contour solution but enables the
+% number of vertices and segments to be preserved.
+% new gB segment vertices are positioned halfway between the original segment
+% midpoint its adjacent segment midpoints.
+%
+%
+% Input
+%  grains - @grain2d - grains output from calcGrains
+%
+% Output
+%  grainsContour - @grain2d with updated grain vertices grains.V
+%
+% Versions
+% 2024-07-23 - Created function, tested on mtexdata twins
+% 2024-07-24 - remove for-loops, tested output is still the same
+
+
+gB = grains.boundary;
+V = gB.V; %initialise V variable as copy of old - this will get updated
+midPoints = gB.midPoint;
+
+x = true(gB.length,1);
+% find gB neighbour segments - attached to either end of the current segment
+% count number of segments adjacent to each segment
+numNeighbours = sum(gB.A_F,2); % A_F is symmetric, doesn't matter whether you sum rows or columns
+% (checked with neighbours1 = sum(gB.A_F,1);)
+
+% if it's not a 'nice' boundary segment with two neighbour segments
+% it's probably a map edge or triple junction or quad point
+% ignore this point
+x(numNeighbours~=2)= false;
+
+% find the positions of these neighbours
+%neihgoburs should end up as a cell array of lenght of num segments
+% neighbours=cell(gB.length,1);
+[n1,nIx,~] = find(gB.A_F); %nix is in sorted order, n1 is not
+neighbours = mat2cell(n1,groupcounts(nIx));
+
+VIds = gB.F(x,:); %vertices of segment X
+vOld = gB.V(VIds); %vertex positions for segment x
+
+%define new vertices as halfway between current and neighbour midpoints
+% but we don't know which vertex the segment starts/ends at
+% so try both (vNew1 vs vNew2)
+vNew1 = (midPoints(permute(cat(2,neighbours{x}),[2 1])) + midPoints(x)) /2 ;
+% reorder the segments
+vNew2 = flip(vNew1,2);
+% decide which way to update V based on minimum distance to the
+% original
+% i.e. don't reassign  the segments backwards
+vDist1 = sum(norm(vNew1-vOld),2); %these are always positive or 0
+vDist2 = sum(norm(vNew2-vOld),2);
+
+% assign new vertices
+vNew = vector3d.nan(size(vNew1)); %initialise vNew
+vNew((vDist1 < vDist2),:) = vNew1((vDist1 < vDist2),:);
+vNew((vDist2 < vDist1),:) = vNew2((vDist2 < vDist1),:);
+if any(vDist2 == vDist1)
+    %unexpected output. assign to a random one and warn user
+    warning('looks like vNew1 and vNew2 are equidistant to vOld? check V reassignment');
+    vNew((vDist2 == vDist1),:) = vNew1((vDist2 == vDist1),:);
+end
+
+% % now update gB.V with vNew
+V(VIds) = vNew; %some things get updated more than once but maybe this is ok
+
+
+grainsContour = grains;
+grainsContour.V = V;

EBSDAnalysis/@grain2d/smooth.m

@@ -1,83 +1,219 @@
-function grains = smooth(grains,iter,varargin)
+function [grains, stablefraction] = smooth(grains,iter,varargin)
 % constraint laplacian smoothing of grain boundaries
 %
 % Input
 %  grains - @grain2d
 %  iter   - number of iterations (default: 1)
 %
+% Optional Input
+%  ebsd   - @EBSD - stop moving gb segments about to cross over to the wrong side of
+% grains.boundary.ebsdId
+%
 % Output
 %  grains - @grain2d
+%  stablefraction - scalar - fraction of boundary vertices that stopped
+%  moving before the last smoothing iteration
 %
 % Options
 %  moveTriplePoints  - do not exclude triple/quadruple points from smoothing
 %  moveOuterBoundary - do not exclude outer boundary from smoothing
 %  second_order, S2  - second order smoothing
-%  rate              - default smoothing kernel  
-%  gauss             - Gaussian smoothing kernel  
-%  exp               - exponential smoothing kernel  
-%  umbrella          - umbrella smoothing kernel   
-
+%  rate              - default smoothing kernel
+%  gauss             - Gaussian smoothing kernel
+%  exp               - exponential smoothing kernel
+%  umbrella          - umbrella smoothing kernel
+
+% grab optional input
+[ebsd,varargin] = getClass(varargin,'EBSD');
+if  ~isempty(ebsd)
+    keepEbsdId=true;
+
+else
+    keepEbsdId=false;
+    stablefraction = [];
+end
+
 if nargin < 2 || isempty(iter), iter = 1; end
 
+if abs(dot(grains.N,zvector)) ~= 1
+    [grains,rot] = rotate2Plane(grains);
+    if keepEbsdId
+        [ebsd] = rotate(ebsd,rot);
+        [grains,stablefraction] = smooth(grains,iter,[{ebsd},varargin]);
+    else
+        grains = smooth(grains,iter,varargin);
+    end
+    grains = inv(rot) * grains;
+    %don't need to backrotate EBSD since it's not returned
+    return
+end
+
+
 % compute incidence matrix vertices - faces
 I_VF = [grains.boundary.I_VF,grains.innerBoundary.I_VF];
 
 % compute vertices adjacency matrix
 A_V = I_VF * I_VF';
 t = size(A_V,1);
+numF = size(I_VF,2); %number of faces (gb segments)
 
 % do not consider triple points
 if check_option(varargin,'moveTriplePoints')
-  ignore = false(size(A_V,1),1);
+    ignore = false(size(A_V,1),1);
 else
-  ignore = full(diag(A_V)) > 2;
+    ignore = full(diag(A_V)) > 2;
 end
 
 % ignore outer boundary
 if ~check_option(varargin,'moveOuterBoundary')
-  ignore(grains.boundary.F(any(grains.boundary.grainId==0,2),:)) = true;
+    ignore(grains.boundary.F(any(grains.boundary.grainId==0,2),:)) = true;
 end
 
 if check_option(varargin,{'second order','second_order','S','S2'})
-  A_V = logical(A_V + A_V*A_V);
-  A_V = A_V - diag(diag(A_V));
+    A_V = logical(A_V + A_V*A_V);
+    A_V = A_V - diag(diag(A_V));
 end
 
 weight = get_flag(varargin,{'gauss','expotential','exp','umbrella','rate'},'rate');
 lambda = get_option(varargin,weight,.5);
 
-V = full(grains.V);
+V = grains.V.xyz;
 isNotZero = ~all(~isfinite(V) | V == 0,2) & ~ignore;
 
+%%%%% this is where the changed bit starts (wrt standard mtex/smooth)
+if keepEbsdId
+    ebsdIdPairs = [grains.boundary.ebsdId; grains.innerBoundary.ebsdId]; % 1 per segment F
+
+    % index me using I_VF -- each column is V belonging to 1 F
+    % so for segment f1
+    % [v1] = find(I_VF(f1,:));
+    %always 2 V per segment F but cellfun doesn't like uniformoutput bc the matrices
+    %are shaped wrong, so reshape stuff afterwards
+    VperF= cellfun(@find,mat2cell(I_VF,size(I_VF,1),ones(1,size(I_VF,2))),'UniformOutput',false); VperF=permute(cat(2,VperF{:}),[2 1]);
+
+    %initialise line segment intersections - this determines whether or not a gb segment may move or not
+    % all segments are allowed to move at the start.
+    % If the the gB segment line intersects the line drawn between the ebsd map
+    % points either size of the gB (ebsdIdPairs), it means that the boundary has
+    % not moved into the wrong part of the ebsd map and you can allow further
+    % smoothing.
+    tfIntersect=true(numF,1); % if tfIntsersect(i) == true, then allow V(i) to move
+end
+
 for l=1:iter
-  if ~strcmpi(weight,'rate')
-    [i,j] = find(A_V);
-    d = sqrt(sum((V(i,:)-V(j,:)).^2,2)); % distance
-    switch weight
-      case 'umbrella'
-        w = 1./(d);
-        w(d==0) = 1;
-      case 'gauss'
-        w = exp(-(d./lambda).^2);
-      case {'expotential','exp'}
-        w = lambda*exp(-lambda*d);
+    if ~strcmpi(weight,'rate')
+        [i,j] = find(A_V);
+        d = sqrt(sum((V(i,:)-V(j,:)).^2,2)); % distance
+        switch weight
+            case 'umbrella'
+                w = 1./(d);
+                w(d==0) = 1;
+            case 'gauss'
+                w = exp(-(d./lambda).^2);
+            case {'expotential','exp'}
+                w = lambda*exp(-lambda*d);
+        end
+
+        A_V = sparse(i,j,w,t,t);
     end
-
-    A_V = sparse(i,j,w,t,t);
-  end
-
-  % take the mean over the neigbours
-  Vt = A_V*V;
-
-  m = sum(A_V,2);
-
-  dV = V(isNotZero,:)-bsxfun(@rdivide,Vt(isNotZero,:),m(isNotZero,:));
-
-  isZero = any(~isfinite(dV),2);
-  dV(isZero,:) = 0;
-
-  V(isNotZero,:) = V(isNotZero,:) - lambda*dV;
-
+
+    % take the mean over the neigbours
+    Vt = A_V * V;
+
+    m = sum(A_V,2);
+
+    dV = V(isNotZero,:)-bsxfun(@rdivide,Vt(isNotZero,:),m(isNotZero,:));
+
+    isZero = any(~isfinite(dV),2);
+    dV(isZero,:) = 0;
+    %%%%%%
+    if keepEbsdId
+        % update V - but save a copy VOld before updating
+        VOld=V;
+    end
+    V(isNotZero,:) = V(isNotZero,:) - lambda*dV;
+    if keepEbsdId
+        %find ebsdIdLines on either side of gB segments F
+        %look through rows for any bad ebsdId points - 0 or other weird numbers
+        badEbsdId = any((~(ebsdIdPairs) | isnan(ebsdIdPairs) |  isinf(ebsdIdPairs)),2);
+        %can't draw line, probably a map edge, exclude me
+        tfIntersect(badEbsdId) = false;
+
+        e1 = ebsd(id2ind(ebsd,ebsdIdPairs(tfIntersect,:))).pos;
+        ebsdIdLine = permute(cat(3,e1.x, e1.y), [2 3 1]);
+
+        %2 Vs per F
+        % test ebsdIdPairs against V
+        % see warning just after the iter for-loop
+        vS= cat(3,V(VperF(tfIntersect,1),1),V(VperF(tfIntersect,1),2)); % xy of segment starts
+        vE= cat(3,V(VperF(tfIntersect,2),1),V(VperF(tfIntersect,2),2)); % xy of segment ends
+        vLine = permute(cat(2,vS,vE),[2 3 1]); %handle dims
+        %if these line segments intersect, this point may still move
+        %if they don't intersect, this V shouldn't move anymore
+        % only repopulate 'true' elements of tfFintersect
+        tfIntersect(tfIntersect) = testIntersection(vLine,ebsdIdLine);
+
+        % if the line segments don't intersect,
+        % revert the vertex positions - don't update V
+        % (just now we did VOld=V; before updating V(isNotZero,:) =
+        % V(isNotZero,:) - lambda*dV;)
+        vDontMove= unique(VperF(~tfIntersect)); % V might appear more than once
+        V(vDontMove,:) = VOld(vDontMove,:);
+
+        % stop moving these vertices in the future
+        % update isNotZero
+        isNotZero(vDontMove) = false;
+    end
+end
+if keepEbsdId
+    % output grain vertices that stopped moving
+    stablefraction = numel(vDontMove)/t;
+end
+% update output
+grains.V = vector3d.byXYZ(V);
+
+end %function
+
+function tf = testIntersection(line1, line2)
+% test whether or not two sets of finite length line segments intersect
+%
+% Inputs
+% line1 = 2*2*n numeric matrix [startx starty; endx endy] * n lines
+% line2 = another line / set of lines similar to line1
+%
+% Outputs
+% tf = 1*n logical array, true if the nth line1 and line2 intersect, false if they don't (or it's another special
+% case*)
+%
+% * if line1==line2 or any of the line lengths are 0
+% (i.e. [startx starty] == [endx endy]), then tf returns
+% false even if they intersect which is ok for us, this
+% should never happen for gb segments which means something has
+% probably gone wrong anyway and we shouldn't move it.
+%
+% References
+% method from here:
+% https://stackoverflow.com/questions/3838329/how-can-i-check-if-two-segments-intersect
+% which references https://bryceboe.com/2006/10/23/line-segment-intersection-algorithm/
+%translation of python code
+% def ccw(A,B,C):
+% return (C.y-A.y) * (B.x-A.x) > (B.y-A.y) * (C.x-A.x)
+%
+% # Return true if line segments AB and CD intersect
+% def intersect(A,B,C,D):
+% return ccw(A,C,D) != ccw(B,C,D) and ccw(A,B,C) != ccw(A,B,D)
+% end
+%
+% def ccw(A,B,C):
+
+aa = permute(line1(1,:,:),[2 3 1]); %startpoint of line1, [x y] * n lines
+bb = permute(line1(2,:,:),[2 3 1]); %endpoint of line1, [x y] * n lines
+cc = permute(line2(1,:,:),[2 3 1]);  %startpoint of line2, [x y] * n lines
+dd = permute(line2(2,:,:),[2 3 1]);  %endpoint of line2, [x y] * n lines
+% use permute to get rid of leading 1* dimension
+
+ccw = @(A,B,C) (C(2,:) - A(2,:)) .* (B(1,:)-A(1,:)) > (B(2,:)-A(2,:)) .* (C(1,:)-A(1,:));
+tf = (ccw(aa,cc,dd) ~= ccw(bb,cc,dd)) & (ccw(aa,bb,cc) ~= ccw(aa,bb,dd));
+tf=tf(:);
 end
 
-grains.V = V;

userScripts/Vivian/smoothTest.m

@@ -0,0 +1,59 @@
+%% WP2 - grain boundary drawing
+% test boundaryContours(grains) and smooth1(grains) functions
+% boundaryContours() redraws grain boundary segments using a contouring 
+% algorithm instead of following pixel edges
+% smooth1() implements constrained laplacian smoothing of grain boundaries
+% with a stopping criteria for based on boundary ebsdId positions
+%
+% Vivian Tong
+% MATLAB R2023a, local MTEX version inside EBSD Interfaces project folder
+%
+% version 2024-07-23 - created and tested on mtexdata twins
+% version 2024-07-24 - remove for-loops in smooth1 and boundaryContours
+
+clear; close all; home
+
+%% Load demo dataset
+ebsd = mtexdata('twins');
+
+%calculate grains
+[grains,ebsd.grainId] = calcGrains(ebsd('indexed'));
+
+%% Describe relevant gb properties
+% check you've understood the gB properties correctly
+%{
+% gB.ebsdInd;  -  EBSD IDs of neighbour segments
+% gB.midPoint; - boundary midpoints
+% gB.V       - [x,y] list of vertices - 3723 long
+% gB.F       - [v1,v2] list of boundary segments (faces)- 3285 long
+%}
+
+%% replace V positions using halfway to midpoints of adjacent Fs
+% similar but not exactly the same as marching squares solution, but keeps
+% existing structure of F and V
+
+tic; grainsContour = boundaryContours(grains); toc;
+
+% % regression test when taking out for-loop
+% tic; grainsContour2 = boundaryContoursFor(grains); toc;
+% figure; plot(ebsd, label2rgb(ebsd.grainId,'jet','w')); hold on;  plot(grainsContour2.boundary,'lineColor','w','linewidth',2);
+% plot(grainsContour.boundary,'lineColor','b'); 
+%% Smooth boundaries
+[grainsContourSmooth10,stablefraction10] = smooth1(grainsContour,10,ebsd('indexed'));
+[grainsContourSmooth30,stablefraction30] = smooth1(grainsContour,30,ebsd('indexed'));
+[grainsContourSmooth100,stablefraction100] = smooth1(grainsContour,100,ebsd('indexed'));
+
+%% Plot stuff
+figure; plot(ebsd, label2rgb(ebsd.grainId,'jet','w','shuffle')); hold on; plot(grains.boundary);
+figure; plot(ebsd, label2rgb(ebsd.grainId,'jet','w','shuffle')); hold on; plot(grainsContour.boundary);
+figure; plot(ebsd, label2rgb(ebsd.grainId,'jet','w','shuffle')); hold on; plot(grainsContourSmooth10.boundary);
+figure; plot(ebsd, label2rgb(ebsd.grainId,'jet','w','shuffle')); hold on; plot(grainsContourSmooth30.boundary);
+figure; plot(ebsd, label2rgb(ebsd.grainId,'jet','w','shuffle')); hold on; plot(grainsContourSmooth100.boundary);
+
+
+figure; histogram(grains.boundary.direction);
+figure; histogram(grainsContour.boundary.direction);
+figure; histogram(grainsContourSmooth10.boundary.direction);
+figure; histogram(grainsContourSmooth30.boundary.direction);
+figure; histogram(grainsContourSmooth100.boundary.direction);
+