diff --git a/EBSDAnalysis/@grain2d/boundaryContours.m b/EBSDAnalysis/@grain2d/boundaryContours.m new file mode 100644 index 000000000..c620ace19 --- /dev/null +++ b/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; diff --git a/EBSDAnalysis/@grain2d/smooth.m b/EBSDAnalysis/@grain2d/smooth.m index cb792e81f..fcb0d3955 100644 --- a/EBSDAnalysis/@grain2d/smooth.m +++ b/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; diff --git a/userScripts/Vivian/smoothTest.m b/userScripts/Vivian/smoothTest.m new file mode 100644 index 000000000..e84bacac5 --- /dev/null +++ b/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); +