diff --git a/EBSDAnalysis/@EBSD/load.m b/EBSDAnalysis/@EBSD/load.m index a7126438a..e569f8c24 100644 --- a/EBSDAnalysis/@EBSD/load.m +++ b/EBSDAnalysis/@EBSD/load.m @@ -95,6 +95,8 @@ end % combine multiple inputs +% ISSUE -- we don't always want to do this e.g. multiple maps in h5 +% file may not be related ebsd = [ebsd{:}]; % ensure unique phases diff --git a/EBSDAnalysis/@grain2d/boundaryContours.m b/EBSDAnalysis/@grain2d/boundaryContours.m new file mode 100644 index 000000000..9578516b4 --- /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.allV +% +% Versions +% 2024-07-23 - Created function, tested on mtexdata twins +% 2024-07-24 - remove for-loops, tested output is still the same +% 2025-03-14 - compatiblity with MTEX 6.2.beta.3 -- switch V = gB.V to gB.allV + +gB = grains.boundary; +V = gB.allV; %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.allV(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.allV with vNew +V(VIds) = vNew; %some things get updated more than once but maybe this is ok + + +grainsContour = grains; +grainsContour.allV = V; diff --git a/EBSDAnalysis/@grain2d/smooth.m b/EBSDAnalysis/@grain2d/smooth.m index 80339b071..18b0956fd 100644 --- a/EBSDAnalysis/@grain2d/smooth.m +++ b/EBSDAnalysis/@grain2d/smooth.m @@ -1,5 +1,9 @@ -function grains = smooth(grains,iter,varargin) -% constraint laplacian smoothing of grain boundaries +function [grains,stablefraction] = smooth(grains,iter,ebsd,varargin) +% constrained laplacian smoothing of grain boundaries +% stopping criteria for smoothing iterations based on boundary ebsdId positions +% +% start excluding gb segments when they cross over to the wrong side of +% grains.boundary.ebsdId % % Input % grains - @grain2d @@ -7,22 +11,33 @@ % % 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 +% Options +% moveTriplePoints - do not exclude triple/quadruple points from +% smoothing (checked - ok with new smooth) +% % the options below have not been checked for compatibility with new smooth! % 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 +% +% +% Versions +% 2024-07-23 - Created function based on mtex/grain2d/smooth +% 2024-07-23 - Replace for-loop with array indexing when excluding segments + if abs(dot(grains.N,zvector)) ~= 1 + % update this to rotate ebsd to plane as well and run smooth1 + [grains,rot] = rotate2Plane(grains); + [ebsd] = rotate(ebsd,rot); + [grains,stablefraction] = smooth1(grains,iter,ebsd,varargin); + grains = inv(rot) * grains; - [grains,rot] = rotate2Plane(grains); - grains = smooth(grains); - grains = inv(rot) * grains; %#ok - return + return end @@ -35,59 +50,165 @@ % compute vertices adjacency matrix A_V = I_VF * I_VF'; t = size(A_V,1); - +numF = size(I_VF,2); % 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 = grains.allV.xyz; -isNotZero = ~all(~isfinite(V) | V == 0,2) & ~ignore; + +V = grains.allV.xyz; % 1 per vertex V +isNotZero = ~all(~isfinite(V) | V == 0,2) & ~ignore; %this is at the start, update as boundary Vs get excluded + +%%%%% this is where the changed bit starts (wrt standard mtex/smooth) +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 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); + + %%%%% from orignal code + 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 neighbors - 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; + %%%%%% + + % update V - but save a copy VOld before updating + VOld=V; + V(isNotZero,:) = V(isNotZero,:) - lambda*dV; + + %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 -grains.allV = vector3d.byXYZ(V,grains.how2plot); +% output grain vertices that stopped moving +stablefraction = numel(vDontMove)/t; + +% update output +grains.allV = vector3d.byXYZ(V); + +end + + +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 diff --git a/interfaces/loadEBSD_oh5.m b/interfaces/loadEBSD_oh5.m new file mode 100644 index 000000000..079d8a5cd --- /dev/null +++ b/interfaces/loadEBSD_oh5.m @@ -0,0 +1,238 @@ +function [ebsd] = loadEBSD_oh5(fname,varargin) +% read EBSD data from EDAX OIM oh5 file +% EBSD data file generated from OIM Analysis software +% can include multiple EBSD maps +% modified from loadEBSD_h5oina +% +% Syntax +% Import EBSD data +% ebsd = loadEBSD_oh5('ebsdData.oh5'); +% +% Flags +% convertSpatial2EulerReferenceFrame +% convertEuler2SpatialReferenceFrame +% +% Options +% +% +% TODO +% 1) option to import only specific EBSD maps +% 2) import other kinds of data (images, EDS?) + +% check file extension before proceeding, avoid mixing up with other hdf5 +% formats +assertExtension(fname,'.oh5'); +if check_option(varargin,'check') + ebsd = []; + return +end + +genericNameFix = {{'\W'},{'_'}}; + +% task: find all groups called EBSD. they may be nested in different +% layers. +locEbsd = locFindH5Groups(fname,'EBSD'); + +if length(locEbsd) > 1 + ebsd = cell(length(locEbsd),1); + disp('more than 1 EBSD dataset in the file, output will be a cell') +end + +for k = 1 :length(locEbsd) % TODO: find a good way to write out multiple datasets + % % EBSD dataset + % EBSD orientations are contained in Group /*/OIM Map*/EBSD/ANG/DATA + EBSD_data = h5info(fname,[locEbsd(k).Name '/Data']); + %EBSD header + EBSD_header = h5info(fname,[locEbsd(k).Name '/Header']); + + % read this EBSD map + % read all EBSD data + EBSDdata = struct; + for thing = 1:length(EBSD_data.Datasets) + sane_name = regexprep(EBSD_data.Datasets(thing).Name,genericNameFix{:}); + EBSDdata.(sane_name)=double(h5read(fname,[EBSD_data.Name '/' EBSD_data.Datasets(thing).Name])); + end + + %read EBSD header + EBSDheader = struct; + for thing = 1:length(EBSD_header.Datasets) + sane_name = regexprep(EBSD_header.Datasets(thing).Name,genericNameFix{:}); + content = h5read(fname,[EBSD_header.Name '/' EBSD_header.Datasets(thing).Name]); + if any(size(content) ~=1) & isnumeric(content) + content = reshape(content,1,[]); + end + EBSDheader.(sane_name) = content; + end + + % find crystal/cartesian alignment from header + crys2cart = split(EBSDheader.Cartesian_Alignment,', '); + xyzAlign = strings(length(crys2cart),3); + xyzAlign(:,1) = extractBefore(crys2cart," || x"); + xyzAlign(:,2) = extractBefore(crys2cart," || y"); + xyzAlign(:,3) = extractBefore(crys2cart," || z"); + % need to adjust the strings so you get a 3*1 string array + % convert to empty strings so that join works + % but reset any missing parts to missing to use later on + xyzAlign(ismissing(xyzAlign))=""; + xyzAlign = join(xyzAlign,"",1); + xyzAlign(xyzAlign=="")=missing; + xyzAlignChars = cellstr(["X||","Y||","Z||"] + xyzAlign); + xyzAlignChars(ismissing(xyzAlignChars)) = []; + + % now set up phases in CSList + phases = h5info(fname,[EBSD_header.Name '/Phase']); + % ---------------- + + CS=cell(1,length(phases.Groups)+1); + CS{1}='notIndexed'; + for phaseN = 1:length(phases.Groups) + pN = ['phase_' num2str(phaseN)]; + EBSDphases.(pN)= struct; + for j = 1:length(phases.Groups(phaseN).Datasets) + sane_name = regexprep(phases.Groups(phaseN).Datasets(j).Name,genericNameFix{:}); + content = h5read(fname,[phases.Groups(phaseN).Name '/' phases.Groups(phaseN).Datasets(j).Name]); + EBSDphases.(pN).(sane_name) = content; + end + + ldims = double([EBSDphases.(pN).Lattice_Constant_a,... + EBSDphases.(pN).Lattice_Constant_b, ... + EBSDphases.(pN).Lattice_Constant_c]); + + % the angle comes in single precision. make sure something + % sufficiently close to 90 resp. 120 does not end up with + % rounding errors instead of using the 'force' option + langle = double(degree*[EBSDphases.(pN).Lattice_Constant_alpha,... + EBSDphases.(pN).Lattice_Constant_beta, ... + EBSDphases.(pN).Lattice_Constant_gamma]); + % PGsymID is not recognised as point group number in MTEX but + % SpaceGroupNumber seems to match SpaceId + % or alternatively Symmetry also matches point group number + csm = crystalSymmetry('SpaceId', EBSDphases.(pN).SpaceGroupNumber); + + if strcmp(csm.lattice,'trigonal') | strcmp(csm.lattice,'hexagonal') + langle(isnull(langle-2/3*pi,1e-7))=2/3*pi; + else + langle(isnull(langle-pi/2,1e-7))=pi/2; + end + + CS{phaseN+1} = crystalSymmetry(csm.pointGroup, ... + ldims,... + langle,... + 'Mineral',char(EBSDphases.(pN).MaterialName),... + xyzAlignChars{:}); + end + + % extract phase names + phaseNames = cellfun(@(x) string(x.mineral),CS(2:end)); + + % fix conflicting phase names + phaseNames = makeDisjoint(phaseNames); + + % write back to CS + for i = 2:length(CS), CS{i}.mineral = char(phaseNames(i-1)); end + + % write out first EBSD dataset + rot = rotation.byEuler(EBSDdata.Phi1,EBSDdata.Phi,EBSDdata.Phi2); + + if size(rot,2)>1 + rot = transpose(rot); + end + + % what data should we read, all or just the standard + phase = EBSDdata.Phase; + pos = vector3d(EBSDdata.X_Position,EBSDdata.Y_Position,0); + + opt=struct; + optList_std = {'IQ' 'CI' 'SEM_Signal' 'Fit' 'PRIAS_Bottom_Strip', 'PRIAS_Center_Square', 'PRIAS_Top_Strip'}; + optNames_std = {'iq', 'ci', 'sem', 'fit', 'PRIAS_Bottom_Strip', 'PRIAS_Center_Square', 'PRIAS_Top_Strip'}; + + + optList = fieldnames(EBSDdata); % all names + %throw out all fields that have already been read in separately + optList = optList(~strcmp('Phase', optList)); + optList = optList(~strcmp('Phi', optList)); + optList = optList(~strcmp('Phi1', optList)); + optList = optList(~strcmp('Phi2', optList)); + optList = optList(~strcmp('X_Position', optList)); + optList = optList(~strcmp('Y_Position', optList)); + + % check if names need replacement + optNames = optList; + for jj = 1:length(optNames) + if any(strcmp(optNames{jj},optList_std)) + optNames{jj} = optNames_std{strcmp(optNames{jj},optList_std)}; + end + end + + % populate opt + for jj = 1:length(optNames) + opt.(optNames{jj}) = EBSDdata.(optList{jj}); + end + + % allow user defined CS, overriding the above + if check_option(varargin,'CS') + CS = get_option(varargin,'CS'); + end + + ebsdtemp = EBSD(pos,rot,phase,CS,opt); + ebsdtemp.opt.Header = EBSDheader; + + %% correct reference frames + % modified from ang loader + + % change reference frame + rotCor = [... + rotation.byAxisAngle(xvector+yvector,180*degree),... % setting 1 + rotation.byAxisAngle(xvector-yvector,180*degree),... % setting 2 + rotation.byAxisAngle(xvector,180*degree),... % setting 3 + rotation.byAxisAngle(yvector,180*degree)]; % setting 4 + + % read the correction setting from h5 EBSD header + corSetting = double(h5read(fname,[EBSD_header.Name '/Coordinate System/ID'])); + + if check_option(varargin,'convertSpatial2EulerReferenceFrame') + flag = 'keepEuler'; + opt = 'convertSpatial2EulerReferenceFrame'; + elseif check_option(varargin,'convertEuler2SpatialReferenceFrame') + flag = 'keepXY'; + opt = 'convertEuler2SpatialReferenceFrame'; + else + if ~check_option(varargin,'wizard') + warning(['.oh5 files use inconsistent conventions for spatial ' ... + 'coordinates and Euler angles. I''m defaulting to ' ... + '''convertEuler2SpatialReferenceFrame'', i.e. ''keepXY'', to correct this']); + end + flag = 'keepXY'; + opt = 'convertEuler2SpatialReferenceFrame'; + end + + ebsdtemp = rotate(ebsdtemp,rotCor(corSetting),flag); + + % set up how2plot + switch flag + case 'keepXY' + ebsdtemp.how2plot=plottingConvention(-zvector,xvector); + case 'keepEuler' + switch corSetting + case 1 + ebsdtemp.how2plot=plottingConvention(zvector,yvector); + case 2 + ebsdtemp.how2plot=plottingConvention(zvector,-yvector); + case 3 + ebsdtemp.how2plot=plottingConvention(zvector,-xvector); + case 4 + ebsdtemp.how2plot=plottingConvention(zvector,xvector); + end + end + + + %build cell if required + if length(locEbsd) > 1 + ebsd{k} = ebsdtemp; + else + ebsd = ebsdtemp; + end + +end + + diff --git a/interfaces/tools/locFindH5Groups.m b/interfaces/tools/locFindH5Groups.m new file mode 100644 index 000000000..dc3bcaabe --- /dev/null +++ b/interfaces/tools/locFindH5Groups.m @@ -0,0 +1,32 @@ +function [ginfo] = locFindH5Groups(fname,pattern) +% modified from loadEBSD_h5 function [ginfo] = locFindEBSDGroups(fname) +% extended to take any string pattern not just 'EBSD' + +info = h5info(fname,'/'); + +ginfo = struct('Name',{},... + 'Groups',{},... + 'Datasets',{},... + 'Datatypes',{},... + 'Links',{},... + 'Attributes',{}); + +ginfo = locGr(fname, info.Groups,ginfo,pattern); +end + +function [ginfo] = locGr(fname,group,ginfo,pattern) + +if ~isempty(group) + + for k=1:numel(group) + attr = group(k).Attributes; + name = group(k).Name; + + if (~isempty(attr) && check_option({attr.Value},pattern)) || endsWith(name,pattern) + ginfo(end+1) = group(k); + end + + [ginfo] = locGr(fname,group(k).Groups,ginfo,pattern); + end +end +end \ No newline at end of file