From 584588ed4282a7eaf2071e49a60499cfe009cb4b Mon Sep 17 00:00:00 2001 From: Brust Date: Mon, 14 Mar 2022 16:41:06 -0400 Subject: [PATCH 1/4] Added CalcGrainSize and modified fine_twin_gbs Added CalcGrainSize function to \auxillary. Computes the ASTM E2627 (area-based) and ASTM E112-13 (lineal intercept) grain size numbers based on only a single phase (trying to stay consistent to your vision) Modified find_twin_gbs function in \auxillary. Changes were minor, added in 'unitCell' to calculation so MTEX doesn't move GBs; return back input ebsd with updated merged parent-twin grain ids (ebsd.grainId); assign parent orientations to grain structure (before, a lot were nans). This orientation assignment is rooted in an assumption (NOTE 2 in script). --- Mart2Aust/auxillary/CalcGrainSize.m | 482 ++++++++++++++++++++++++++++ Mart2Aust/auxillary/find_twin_gbs.m | 86 +++-- 2 files changed, 536 insertions(+), 32 deletions(-) create mode 100644 Mart2Aust/auxillary/CalcGrainSize.m diff --git a/Mart2Aust/auxillary/CalcGrainSize.m b/Mart2Aust/auxillary/CalcGrainSize.m new file mode 100644 index 0000000..aa5f5b4 --- /dev/null +++ b/Mart2Aust/auxillary/CalcGrainSize.m @@ -0,0 +1,482 @@ +function [varargout] = CalcGrainSize(varargin) + + % Create Grain (and Twin--if necessary) structure for data output + E2627 = []; + Twins = []; + + % First object is always the ebsd + ebsd = varargin{1}; + % Second object is always the phase to extract + PhaseId = varargin{2}; + % Third object is always the flag of whether to merge the twins into + % their parent orientation + MergeTwins = varargin{3}; + + % Fourth object would be the grain boundary misorientation tolerance + % angle (what classifies the GB interface: default is 5 degrees) + if length(varargin) == 4 + GB_Tol = varargin{4} * degree; + else + GB_Tol = 5*degree; + end + + PxTol = 100; + % NOTE 1: Pixel tolerance can also be input argument, although default + % ASTM standards is 100 pixels; + + % Determine number of phases present in ebsd dataset. + EbPh = ebsd(ismember(ebsd.phaseId,PhaseId)); + + % Compute grain boundaries and characteristics for Phase and + % extract all grain ids for Phase + if MergeTwins + [GrnPh,EbPh,S3,S9] = find_twin_gbs(EbPh, PhaseId); + Twins.S3 = S3; + Twins.S9 = S9; + else + [GrnPh,EbPh.grainId] = calcGrains(EbPh,'unitCell','angle',GB_Tol); + end + GrnId = GrnPh.id; + + % Initiate indexing for which grains fail + GrnFail = (GrnPh.grainSize < PxTol); + + % Ebsd step size (ASTM scaling factors) + stepsz = EbPh.scanUnit; + if stepsz == 'um' + delta = 1e-6; + scale = 1e6; + elseif stepsz == 'mm' + delta = 1e-3; + scale = 1; + else + warning('Ebsd step size not recorded') + end + + % Extract the unique grain IDs bordering the IPF map and ignore + % from calculations (along with smaller grains) + GrnBB = GrnPh.isBoundary; + GrnFail(GrnBB) = 1; + + % Also extract grains included in ASTM + GrnPass = GrnId(logical(abs(GrnFail-1))); + + % Extract area and % multiply measured grain area by (step size)^2 + ASTM_GrnArea = area(GrnPh(GrnPass)) * delta^2; + + % If hex-grid, modify to follow ASTM E2627 standards + if length(EbPh.unitCell) == 6 + ASTM_GrnArea = ASTM_GrnArea * sqrt(3)/2; + end + + % Compute ASTM E2627 grain size number + ASTM_E2627_G = -3.3223*log10(mean(ASTM_GrnArea*scale))-2.995; + + % Collect statistics from grain size calculation + Sz = size(GrnPh(GrnPass),1); + freqGS = 1/Sz * sum(ASTM_GrnArea); + StDev = sqrt(1/(Sz-1) * sum((ASTM_GrnArea - freqGS).^2)); + StError = StDev / sqrt(Sz); + CI95 = tstat3(max(Sz-1,1), 0.975, 'inv'); + yCi95 = bsxfun(@times, StError, CI95(:)); + + % Store ASTM E2627 Grain Size Calculations + E2627.AcceptableGrainIds = GrnPass; + E2627.GrainSize = ASTM_E2627_G; + E2627.Ci95 = yCi95; + + % Calculate the ASTM E112_13 grain size number + [E112_13] = Calculate_ASTM_E112_13(GrnPh,EbPh); + + % Pack the kids up and go home + varargout{1} = EbPh; + varargout{2} = GrnPh; + varargout{3} = E2627; + varargout{4} = E112_13; + if isempty(Twins) == 0 + varargout{5} = Twins; + end + +end + +function [E112_13] = Calculate_ASTM_E112_13(varargin) + +%% Function Description + +% This function follows ASTM E112-13 standards (for grains) by combining +% annealling twins with the parent s.t. the parent-twin system becomes a +% single entity whose grain boundaries are defined by the parent grain. + +% Single intercept lines are performed on any packet, block, and sub-block +% boundaries for further analysis if the material is steel by randomly +% assigning a start and stop on the x_min and y_max edge, and then +% performing the same sort of standard segmentation analysis. + +% The grain boundary lineal intercept method involves 12 total lineal +% intercepts, with 8 vertical lines spanning the x-axis coordinates and +% 4 horizontal lines spanning the y-axis, all lines being equally +% distributed from each other. + + +%% Function Implementation + + % Grab grains and ebsd + grains = varargin{1}; + ebsd = varargin{2}; + + % Assign myEBSD and Grain structures (called Feature here) + Bounds = grains.boundary; + + % If user wants to define their own line segment, use their + % coordinates. Else, randomly confine a line to the minimum x value + % (left edge), a y value contained within this edge, and two randomly + % chosen x,y endpoints. + if length(varargin)==3 + X1 = varargin{3}(1); + X2 = varargin{3}(2); + Y1 = varargin{3}(3); + Y2 = varargin{3}(4); + + % Establish the line endpoints + L1 = [X1,Y1]; + L2 = [X2,Y2]; + % Length of line + lenLine = sqrt((L2(1)-L1(1))^2+(L2(2)-L1(2))^2); + + % Determine number of intersections and line-segment lengths + [X,Y,SortedLineSegs] = Bounds.intersect(L1,L2); + Xint = X(find(~isnan(X))); + Yint = Y(find(~isnan(Y))); + % Number of intersections + P_i = length(SortedLineSegs); + + % Average grain size + li = lenLine / P_i; + % Surface Area + Sv = 2/lenLine*1000; + + % Calculate ASTM Grain Size Number Per ASTM E112-13 + ASTM_Num = -6.643856*log10(lenLine/1000)-3.288; + + else % Do multiple lineal intercepts for ASTM grain size measurement + + % Extract triple point coordinates and set pixelated tolerance + tpX = Bounds.triplePoints.x; + tpY = Bounds.triplePoints.y; + tpTol = 1; + + % Range of x and y values and lineal intercept line length + Lx = range(ebsd.x); + Ly = range(ebsd.y); + TotLine = (Lx*4)+(Ly*8); + + % X-based values + XPart = Lx / 10; + XStrt = min(ebsd.x); + XEnd = max(ebsd.x); + + % Y-based values + YPart = Ly / 6; + YStrt = min(ebsd.y); + YEnd = max(ebsd.y); + + % X-dependent line that runs across the x-axis coordinates at + % constant y values + Xdep_X1 = ones(1,4)*XStrt-10; + Xdep_X2 = ones(1,4)*XEnd+10; + Xdep_Y1 = linspace(YStrt+YPart,YEnd-YPart,4); + Xdep_Y2 = Xdep_Y1; + + % Y-dependent line that runs across the y-axis coordinates at + % constant x values + Ydep_X1 = linspace(XStrt+XPart,XEnd-XPart,8); + Ydep_X2 = Ydep_X1; + Ydep_Y1 = ones(1,8)*YStrt-10; + Ydep_Y2 = ones(1,8)*YEnd+10; + + % Number of lineal intercepts + Sz = 12; + + % Preallocate arrays + li = zeros(Sz,1); % Line segment length + NSegs = zeros(Sz,1); % Number of intersections + Sv = zeros(Sz,1); % Surface Area + ASTM_Nums = zeros(Sz,1); % ASTM E112-13 grain size number + + % Loop through all line segments and calculate grain size based on + % lineal intercept method + for ii = 1:Sz + + if ii > 4 + X1 = Ydep_X1(ii-4); + X2 = Ydep_X2(ii-4); + Y1 = Ydep_Y1(ii-4); + Y2 = Ydep_Y2(ii-4); + else + X1 = Xdep_X1(ii); + X2 = Xdep_X2(ii); + Y1 = Xdep_Y1(ii); + Y2 = Xdep_Y2(ii); + end + + % Establish the line endpoints + L1 = [X1,Y1]; + L2 = [X2,Y2]; + + % Determine number of intersections and line-segment lengths + [X,Y,~] = Bounds.intersect(L2,L1); + + % Extract the intersections + Ints = ~isnan(X); + % Id grains of the intersections. + IntGrns = Bounds.grainId(Ints,:); + % NOTE 2: The first and last entries will always be 0 and an + % existing grain id. These first 0's are the edges of the + % boundaries and count for 0.5 intersections + + % Determine first whether any repeated grain intersection Ids + % coincide with the ends of the micrograph + RepeatGrnEnds = zeros(length(IntGrns),1,'logical'); + EndIds = [1,length(IntGrns)]; + % Find whether any "middle" grain intersection pattern repeats + % one of the ends of the micrograph ids. + MiddleGrnRepeat = ismember(IntGrns(2:EndIds(2)-1,:),IntGrns(EndIds,:),'row'); + RepeatGrnEnds(2:EndIds(2)-1) = MiddleGrnRepeat; + % NOTE 3: If this does occur, it means that either an unindexed + % region or different phase grain is sandwiched between one of + % end grains and another grain. We need to delete these grains + % and the corresponding line segments for accurate line lengths + + % Delete those grains that repeat at the end + IntGrns(RepeatGrnEnds,:) = []; + + % Find X and Y coordinates of intersections + Xint = X(Ints); + Yint = Y(Ints); + + % Extract coordinates of the unique grains at intersections + [UnqIntGrns,UnqIntGrnId,UnqIntRepIds] = unique(IntGrns,'rows','stable'); + UnqXint = Xint(UnqIntGrnId); + UnqYint = Yint(UnqIntGrnId); + + % Look for tangents (where line runs across 2 grain boundaries + % and intersects multiple times) + RepCnts = histcounts(UnqIntRepIds); + RepGrns = UnqIntGrns(RepCnts > 2,:); + RepGrnsRev = fliplr(RepGrns); + % If the reverse case also exists [i.e. (i,j) and (j,i), + % # delete (j,i) from the unique grain intersection list + if any(ismember(RepGrnsRev,UnqIntGrns,'rows')) + DelGrnRow = ismember(UnqIntGrns,RepGrnsRev,'rows'); + UnqIntGrns(DelGrnRow,:) = []; + end + + % Determine the line segment lengths. Pick from the second + % column because it ignores all cases where a grain intersects + % with an unindexed or different phase grain + UnqGrnSegs = unique(UnqIntGrns(:,2),'stable'); + lenline = 0; + for IntCnt = 1:length(UnqGrnSegs) + [~,~,lenseg] = grains(UnqGrnSegs(IntCnt)).boundary.intersect(L2,L1); + if length(lenseg) > 1 + lenseg = sum(lenseg); + end + lenline = lenline+lenseg; + end + + % Those flagged as tangents will now reduce the intersection + % count by 0.5x tangent grain boundaries + TanFlgs = ismember(UnqGrnSegs,RepGrns); + % NOTE 4: ASTM counts tangent boundaries as 1 intersection, + % whereas MTEX will flag as 2, so we flag here and then + % subtract 0.5 intersections for eacht time this occurs. Ergo, + % grain i and j show a tangent boundary, initially this would + % count for 2 intersections, but each grain will get flagged + % and thus 2-0.5*2 = 1 intersection. + % NOTE 5: ASTM doesn't specify, but here, if a grain is tangent + % with an non-phase grain boundary, that intersection only + % counts as 0.5 total + + % Sort based on X intercepts unless x is constant, then sort by Y + if X2 - X1 == 0 + [~,SortIntGrnId] = sort(UnqYint,'descend'); + else + [~,SortIntGrnId] = sort(UnqXint); + end + + % Extract the sorted, unique grain ids of intercepts and x and y + % coordinate intercepts. Although sorting is unneccessary, it + % allows us to follow the line from grain to grain either + % horizontally or vertically + SortedXint = UnqXint(SortIntGrnId)'; + SortedYint = UnqYint(SortIntGrnId)'; + + % Number of intersections. A full intersection consists of a + % line segment passing through two same-phase grains; if a 0 + % exists, it's either a boundary grain or a different phase + % (unindexed). These intersections count as 0.5 intersections. + P_i = length(UnqIntGrns) / 2 - 0.5*sum(TanFlgs); + + % Check for triple points + tpDsts = abs((L2(2)-L1(2)).*tpX-(L2(1)-L1(2)).*tpY+L2(1)*L1(2)-L2(2)*L1(1))/... + sqrt((L2(2)-L1(2))^2+(L2(1)-L1(1))^2); + + % For square pixels, it's simple + if length(ebsd.unitCell) == 4 + % Calculate the tolerance for triple point intersection + % distances + tolDst = 2*tpTol*abs(ebsd.unitCell(1,1)); + else + % For hexagonal pixels, keep it simple and only use the + % x-distance between pixel centroids + XminHex = min(ebsd.unitCell(:,1)); + XmaxHex = max(ebsd.unitCell(:,1)); + tolDst = tpTol*(XmaxHex-XminHex); + end + + % Find the length of the number of triple points whose + % tolerance is intersected by the lineal intercept line + tpSegs = length(find(tpDsts < tolDst)); + + % Adjust the number of intersections based on the number of + % triple point intersections + P_i = P_i + tpSegs*0.5; + + % Now check for any lineal intercepts that directly intersect a + % triple point. + TPtChk = zeros(length(SortedXint),1); + % NOTE 6: The way MTEX searches for grain boundaries, triple + % points don't necessarily (nor commonly) land on the vertices + % of indexed pixels, but rather at points within the pixels. + % NOTE 7: Because of this, if a triple point does happen to + % land on the vertex of a pixel, it results in 3 detected + % grain boundary segmentations. Since ASTM standards call for + % 1.5 segmentations at triple points, we need to adjust P_i by + % subtracting 2 (since we already added the 0.5 in the previous + % line) + + % Loop through all of the intersected coordinates and find + % those that happen to align with a triple point coordinate + for kk = 1:length(SortedXint) + TPtChk(kk) = any(tpX == SortedXint(kk) & tpY == SortedYint(kk)); + + end + + % If any triple points exist, modify our intersection count + % by reducing by 2xn, where n is the number of triple points + % intersected + if length(find(TPtChk)) > 0 + P_i = P_i - length(find(TPtChk == 1)); + end + + % Mu line segment length + if length(P_i) > 1 + keyboard + end + li(ii) = lenline / P_i; + % Number of intersections + NSegs(ii) = P_i; + % Surface Area + Sv(ii) = 2/li(ii)*1000; + + % ASTM E112-13 grain size number + ASTM_Nums(ii) = -6.643856*log10(li(ii)/1000)-3.288; + end + + % CI Functions + CI95 = tstat3(max(Sz-1,1), 0.975, 'inv'); + + % Lineal Intercept Data + E112_13.LinealIntercept.Mu = TotLine / sum(NSegs); + E112_13.LinealIntercept.LineLength = TotLine; + E112_13.LinealIntercept.Sigma = std(li); + E112_13.LinealIntercept.StErr = std(li) / sqrt(Sz); + E112_13.LinealIntercept.CI95 =... + bsxfun(@times, std(li)/sqrt(Sz), CI95(:)); + + % ASTM Grain Size Data + E112_13.GrainSize.Nums = ASTM_Nums; + E112_13.GrainSize.Mu = mean(ASTM_Nums); + E112_13.GrainSize.StDev = std(ASTM_Nums); + E112_13.GrainSize.ASTM_StError = std(ASTM_Nums) / sqrt(Sz); + E112_13.GrainSize.ASTM_CI95 =... + bsxfun(@times, std(ASTM_Nums)/sqrt(Sz), CI95(:)); + + % Surface Volume Data + E112_13.SurfaceVolume.Mu = mean(Sv); + E112_13.SurfaceVolume.Sigma = std(Sv); + E112_13.SurfaceVolume.StError = std(ASTM_Nums) / sqrt(Sz); + E112_13.SurfaceVolume.CI95 =... + bsxfun(@times, std(ASTM_Nums)/sqrt(Sz), CI95(:)); + + % Save the number of segmentations + E112_13.NumberOfSegmentations = NSegs; + + end +end + +function pt = tstat3( v, tp, stat ) + % TSTAT3 computes one of three t_statistics: one-sided t-probability, or + % two-sided t-probability, or the inverse t-statistic in any single + % call. It does not take vectors as arguments. + % + % INPUT ARGUMENTS: + % v Degrees-of freedom (integer) + % + % tp T-statistic: to produce t-probability + % OR + % Probability (alpha): to produce t-statistic + % + % stat Desired test ? + % 'one' One-tailed t-probability + % 'two' Two-tailed t-probability + % 'inv' Inverse t-test + % + % OUTPUT ARGUMENT: + % pt T-probability OR T-statistic, as requested in ?stat? + % + % USE: + % FIND ONE-TAILED PROBABILITY GIVEN t-STATISTIC & DEGREES-OF-FREEDOM + % p = tstat3(v, t, 'one') + % + % FIND TWO-TAILED PROBABILITY GIVEN t-STATISTIC & DEGREES-OF-FREEDOM + % p = tstat3(v, t, 'two') + % + % FIND ONE-TAILED t-STATISTIC GIVEN PROBABILITY & DEGREES-OF-FREEDOM + % t = tstat3(v, p, 'inv') + % + % + % + % Star Strider ? 2016 01 24 ? + % % % % T-DISTRIBUTIONS ? + % Variables: + % t: t-statistic + % v: degrees of freedom + tdist2T = @(t,v) (1-betainc(v/(v+t^2), v/2, 0.5)); % 2-tailed t-distribution + tdist1T = @(t,v) 1-(1-tdist2T(t,v))/2; % 1-tailed t-distribution + % This calculates the inverse t-distribution (parameters given the + % probability ?alpha? and degrees of freedom ?v?: + t_inv = @(alpha,v) fzero(@(tval) (max(alpha,(1-alpha)) - tdist1T(tval,v)), 5); % T-Statistic Given Probability ?alpha? & Degrees-Of-Freedom ?v? + statcell = {'one' 'two' 'inv'}; % Available Options + nc = cellfun(@(x)~isempty(x), regexp(statcell, stat)); % Logical Match Array + n = find(nc); % Convert ?nc? To Integer + if (length(v) > 1) || (length(tp) > 1) + error(' ?> TSTAT3 does not take vectorised inputs.') + elseif isempty(n) % Error Check ?if? Block + error(' ?> The third argument must be either ''one'', ''two'', or ''inv''.') + elseif (n == 3) && ((tp < 0) || (tp > 1)) + error(' ?> The probability for ''inv'' must be between 0 and 1.') + elseif (isempty(v) || (v <= 0)) + error(' ?> The degrees-of-freedom (''v'') must be > 0.') + end + switch n % Calculate Requested Statistics + case 1 + pt = tdist1T(tp, v); + case 2 + pt = tdist2T(tp, v); + case 3 + pt = t_inv(tp, v); + otherwise + pt = NaN; + end +end diff --git a/Mart2Aust/auxillary/find_twin_gbs.m b/Mart2Aust/auxillary/find_twin_gbs.m index cd70fb5..875f814 100644 --- a/Mart2Aust/auxillary/find_twin_gbs.m +++ b/Mart2Aust/auxillary/find_twin_gbs.m @@ -1,34 +1,56 @@ -function [Aus_gb_merged,S3,S9] = find_twin_gbs(ebsd, phaseID) -% Returns all the grain boundaries, as well as the S3 and S9 boundaries - -% find all boundaries -[grains,~,~] = calcGrains(ebsd, 'angle', 1*degree); -grains = grains.smooth(5); -Recon_phasename = ebsd.CSList{3}.mineral; -Aus_gb = grains.boundary(Recon_phasename,Recon_phasename); - -% define S3 and S9. -%S3 is defined using u1v1 method from Steve -CS_A = ebsd.CSList{3}; -u1 = Miller( 1, 2, 1,CS_A); -v1 = Miller(-1,-2,-1,CS_A); -u2 = Miller(-1, 0, 1,CS_A); -v2 = Miller( 1, 0,-1,CS_A); -S3_twin_miso = orientation('map',u1,v1,u2,v2); -% S9 was done from 1984 Hans Grimmer Paper -S9_angle = 38.94*degree; -S9_axis = vector3d(1,1,0); -S9_twin_miso = orientation('axis',S9_axis,'angle',S9_angle,CS_A); -S9_twin_miso.SS = S9_twin_miso.CS; - -% restrict to twinnings with threshold 0.5 degree -S3_isTwinning = angle(Aus_gb.misorientation,S3_twin_miso) < 1.5*degree; -S9_isTwinning = angle(Aus_gb.misorientation,S9_twin_miso) < 1.5*degree; - -S3 = Aus_gb(S3_isTwinning); -S9 = Aus_gb(S9_isTwinning); -All_Twin_Boundaries = [S9 S3]; - -[Aus_gb_merged,~] = merge(grains,All_Twin_Boundaries); +function [Aus_gb_merged,ebsd,S3,S9] = find_twin_gbs(ebsd, phaseID) + % Returns all the grain boundaries, as well as the S3 and S9 boundaries, + % and stores the grain id in the ebsd data set + + % find all boundaries + [grains,ebsd.grainId,~] = calcGrains(ebsd,'unitCell','angle', 1*degree); + grains = grains.smooth(5); + Recon_phasename = ebsd.CSList{3}.mineral; + Aus_gb = grains.boundary(Recon_phasename,Recon_phasename); + % NOTE 1: Alex add in: 'unitCell' doesn't allow for grains to + % eat up other grains (phases, unindexed points, etc) + + % define S3 and S9. + %S3 is defined using u1v1 method from Steve + CS_A = ebsd.CSList{3}; + u1 = Miller( 1, 2, 1,CS_A); + v1 = Miller(-1,-2,-1,CS_A); + u2 = Miller(-1, 0, 1,CS_A); + v2 = Miller( 1, 0,-1,CS_A); + S3_twin_miso = orientation('map',u1,v1,u2,v2); + % S9 was done from 1984 Hans Grimmer Paper + S9_angle = 38.94*degree; + S9_axis = vector3d(1,1,0); + S9_twin_miso = orientation('axis',S9_axis,'angle',S9_angle,CS_A); + S9_twin_miso.SS = S9_twin_miso.CS; + + % restrict to twinnings with threshold 0.5 degree + S3_isTwinning = angle(Aus_gb.misorientation,S3_twin_miso) < 1.5*degree; + S9_isTwinning = angle(Aus_gb.misorientation,S9_twin_miso) < 1.5*degree; + + S3 = Aus_gb(S3_isTwinning); + S9 = Aus_gb(S9_isTwinning); + All_Twin_Boundaries = [S9 S3]; + + % Merge parents and twins + [Aus_gb_merged,ParentId] = merge(grains,All_Twin_Boundaries); + + % Loop through ebsd and replace grainIds with merged ones + for ii = 1:length(grains.id) + ebsd.grainId(ismember(ebsd.grainId,ii)) = ParentId(ii); + end + + % Assign parent orientations to the merged grain structure. + for jj = 1:length(Aus_gb_merged.id) + MergedGrns = grains.id(ismember(ParentId,jj)); + [~,MaxId] = max(grains.grainSize(MergedGrns)); + ParOr = grains(MergedGrns(MaxId)).meanOrientation; + Aus_gb_merged(jj).meanOrientation = ParOr; + end + % NOTE 2: When you merged the grains, a lot of the grain orientations + % were assigned values of nan. This fixes that, but under the + % assumption that the largest surface area for each identified parent + % and twin is the parent, and all others are the twins. This isn't + % necessarily correct, but can be fixed later on and is fast end \ No newline at end of file From c3d407bb416bee91667b84e955debbf90d1272a2 Mon Sep 17 00:00:00 2001 From: Brust Date: Mon, 14 Mar 2022 17:14:15 -0400 Subject: [PATCH 2/4] Truncate_Ebsd_Function Lazy ebsd truncation implementation to make large ebsd data sets smaller based on input min,max(x,y) coordinates --- Mart2Aust/auxillary/truncate_ebsd.m | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 Mart2Aust/auxillary/truncate_ebsd.m diff --git a/Mart2Aust/auxillary/truncate_ebsd.m b/Mart2Aust/auxillary/truncate_ebsd.m new file mode 100644 index 0000000..518b023 --- /dev/null +++ b/Mart2Aust/auxillary/truncate_ebsd.m @@ -0,0 +1,7 @@ +function [trunc_ebsd] = truncate_ebsd(ebsd,xmin,xmax,ymin,ymax) + + rr = [xmin ymin (xmax-xmin) (ymax-ymin)]; + cond = inpolygon(ebsd,rr); + trunc_ebsd = ebsd(cond); + +end \ No newline at end of file From 6258c163dcdf4099cd3fe4077fba7d665c75c192 Mon Sep 17 00:00:00 2001 From: Brust Date: Tue, 21 Jun 2022 14:19:04 -0400 Subject: [PATCH 3/4] Updated Functions Updated functions HT_grain_guess (indexing issue when calling calcGrains -- we can discuss further, stupid MTEX nonsense), and then necessary minor changes to functions AutoOR_estimation and Call_Reconstruction --- Mart2Aust/Core/AutoOR_estimation.m | 8 +- Mart2Aust/Core/Call_Reconstruction.m | 33 ++- Mart2Aust/auxillary/HT_grain_guess.m | 310 +++++++++++++-------------- 3 files changed, 184 insertions(+), 167 deletions(-) diff --git a/Mart2Aust/Core/AutoOR_estimation.m b/Mart2Aust/Core/AutoOR_estimation.m index b0504fc..a6ed830 100644 --- a/Mart2Aust/Core/AutoOR_estimation.m +++ b/Mart2Aust/Core/AutoOR_estimation.m @@ -34,7 +34,7 @@ %% search the low temp phase (Phase 2) for acceptably close PAG % for first iteration, make the entire scan searchable -searchable_ebsd = ebsd(ebsd.phaseId == 2); +searchable_ebsd = ebsd; for iteration = 1:10 % NOTE TO FUTURE USERS: The first step in this code is to use the MTEX % calcgrain tool to find likely martensite grains, then merge grains @@ -66,7 +66,7 @@ %to the leftovers-only search, JIC. Feel free to try either way. %% Method 1: Search all unused areas % find indices of a possible grain - [Aus_Grain_Ids,~,AusOr] = HT_grain_guess(searchable_ebsd); + [Aus_Grain_Ids,~,AusOr] = HT_grain_guess_Alex(searchable_ebsd); % extract the expected grain and leave the remainder for future searches grain_ebsd = searchable_ebsd(Aus_Grain_Ids); searchable_ebsd(Aus_Grain_Ids) = []; @@ -89,6 +89,7 @@ martensite=martensite(keep(1:options.OR_sampling_size)); else end + % plot the martensite pole figure and grain if desired if options.OR_plot_PAG_Mart_Guess == 1 figure(); @@ -142,7 +143,6 @@ %leftover can ignore austenite_prior_odf=uniformODF(CS_HT,SS); - prior_pars=struct; prior_pars.ksi_prior_mu=ksi_prior_mu; prior_pars.ksi_prior_sigma=ksi_prior_sigma; @@ -171,7 +171,7 @@ count = count+1; % Optimization function outside of ML add-on [MAPpars,~,~,~]=fminsearch(optimfunc,x0,MAP_options); - + keyboard if (MAPpars(1) < MAPpars(2) && MAPpars(3)) || count > 2 init_guess = 1; % Enforce constraint if need be diff --git a/Mart2Aust/Core/Call_Reconstruction.m b/Mart2Aust/Core/Call_Reconstruction.m index 9e51e54..d25f19c 100644 --- a/Mart2Aust/Core/Call_Reconstruction.m +++ b/Mart2Aust/Core/Call_Reconstruction.m @@ -1,4 +1,4 @@ -function HT_ebsd = Call_Reconstruction(orig_ebsd,options,LT_MDF) +function [HT_ebsd, GlobOP_wts] = Call_Reconstruction(orig_ebsd,options,LT_MDF) % Do the actual reconstruction from the low temp to the high temp phase % Pre-flight stuff: @@ -23,7 +23,7 @@ % Now do the actual reconstruction as a seperate function (helps avoid % accidental overwrites and removes pre-flight parameters) -recon_ebsd = Reconstruct_HT_grains(LT_ebsd,... +[recon_ebsd, RecOP_wts] = Reconstruct_HT_grains(LT_ebsd,... neigh_list,... IP_wts,... options); @@ -34,10 +34,15 @@ HT_ebsd = orig_ebsd; HT_ebsd(orig_ebsd.phaseId == 2).orientations = temp_ebsd.orientations; HT_ebsd(orig_ebsd.phaseId == 2).phaseId = recon_ebsd.phaseId; + +% Set likelihood for entire ebsd data set, regardless of the phase. Those +% specific likelihood values can be later extracted using ebsd('phase').id +GlobOP_wts = zeros(length(orig_ebsd),1); +GlobOP_wts(orig_ebsd.isIndexed) = RecOP_wts; end -function LT_ebsd = Reconstruct_HT_grains(LT_ebsd,neigh_list,IP_wts,options) +function [LT_ebsd,RecOP_wts] = Reconstruct_HT_grains(LT_ebsd,neigh_list,IP_wts,options) %Reconstruct_HT_grains perform the actual reconstruction % performs a 'while ' loop that progressively cuts out possible grains % until eithr there is nothing left to cut, or the algorithm starts @@ -51,6 +56,9 @@ %%figure() %%plot(Active_Ebsd,Active_Ebsd.orientations) +% Alex add on: Global Likelihood +RecOP_wts = zeros(length(LT_ebsd),1); + while continue_recon == 1 % checks to break out of while loop if iterations > options.max_recon_attempts || bad_cut_counter > 10 @@ -61,7 +69,7 @@ disp('=========================================\n') continue end - if sum(LT_ebsd.phaseId == 2) == 0 + if sum(LT_ebsd.phaseId == 2) < 4 continue_recon = false; disp('\n=========================================') disp('reconstruction completed successfully!') @@ -71,7 +79,12 @@ %% ======== Step 1 ======== %% % Find a likely high temp grain orientation to attempt to cut out. - Guess_ID = randsample(1:length(Active_Ebsd),1,true,Active_Ebsd.ci); + try + Guess_ID = randsample(1:length(Active_Ebsd),1,true,Active_Ebsd.ci); + catch + Guess_ID = randsample(1:length(Active_Ebsd),1,true); + end + HT_guess_ori = Active_Ebsd(Guess_ID).orientations; % Using that guess, make a rough cut to find a likely HT grain.If the @@ -102,7 +115,7 @@ % if the grain made it this far, reset the counters and do 5 graph % cuts (1 for the parent, 4 for the twins) to find the parent/twin areas bad_cut_counter = 0; - [PT_ID_stack] = Seperate_Twins(proposed_grain, neigh_list, IP_wts, PT_oris,options); + [PT_ID_stack,LocOP_wts] = Seperate_Twins(proposed_grain, neigh_list, IP_wts, PT_oris,options); % Use this data to overwrite the orientation and phase data in LT_ebsd for i = 1:size(PT_oris,1) mask = PT_ID_stack(i,1:end); @@ -113,6 +126,7 @@ IDs_to_assign = proposed_grain(mask).id; LT_ebsd(ismember(LT_ebsd.id, IDs_to_assign)).orientations = PT_oris(i); LT_ebsd(ismember(LT_ebsd.id, IDs_to_assign)).phase = 3; + RecOP_wts(ismember(LT_ebsd.id, IDs_to_assign)) = LocOP_wts(mask); end % Prune out any orphaned pixels @@ -292,6 +306,7 @@ [~,~,cs,~]=maxflow(FP_digraph,N+1,N+2); cs(cs>length(Rough_Guess)) = []; proposed_grain = Rough_Guess(cs); + % Code for debugging to show the cut out area. NOTE: this is not a grain, % Its just a region of the scan that likely has at least one complete grain % in it @@ -307,7 +322,7 @@ end -function [Parent_and_Twin_IDs] = Seperate_Twins(proposed_grain, neigh_list, IP_wts, PT_oris,options) +function [Parent_and_Twin_IDs,LocOP_wts] = Seperate_Twins(proposed_grain, neigh_list, IP_wts, PT_oris,options) % Given we have a for sure parent grain, check to see if parts would make % more sense as a twin or as part of the parent. @@ -338,6 +353,7 @@ mx = (bg_likelyhoods*options.RGC_post_pre_m); b = options.RGC_post_pre_b; OP_Parent_wts = mx + b; +LocOP_wts = OP_Parent_wts; clear Parent_variants Parent_odf bg_likelyhoods mx b for i = 2:size(PT_oris,1) @@ -375,6 +391,8 @@ % save out the mask of the cut % yes_mask = [1:N].*ismember([1:N],cs); Parent_and_Twin_IDs(i,1:end) = ismember((1:N),cs); + % Establish OP weights for completed parent and twins + LocOP_wts(ismember((1:N),cs)) = OP_twin_wts(cs); % find the already assigned orientations, set their in-plane weights to % zero so they won't get pulled out again. assigned =(1:N).*(sum(Parent_and_Twin_IDs,1)>0); @@ -503,4 +521,3 @@ % neighborhood) change options.degree_of_connections_for_neighborhood % from 2 to 1. - diff --git a/Mart2Aust/auxillary/HT_grain_guess.m b/Mart2Aust/auxillary/HT_grain_guess.m index 89a686a..07ed473 100644 --- a/Mart2Aust/auxillary/HT_grain_guess.m +++ b/Mart2Aust/auxillary/HT_grain_guess.m @@ -2,161 +2,161 @@ % combination of MTEX-driven and graph cut `reconstruction' algorithm used % to segment a portion of a grain expected to be from a single Austenite % grain, which can then be used for the automatic determination of the OR - -ksiKS = [5.26,10.30,10.53]; -ksiNW = [0,9.74,9.74]; -CS_LT = ebsd.CSList{2}; -CS_HT = ebsd.CSList{1}; - -% Allow the use of an initial Orientation Relationship guess. Otherwise, -% default to GT (ask eric, not sure what this is -if exist('OR_guess','var') - assert(isvector(OR_guess) && length(OR_guess) ==3 && isfloat(OR_guess),... - "OR_guess must be a length 3 array of floats") - ksi = OR_guess; % Begin with provided guess -else - ksi = [2.16,8.06,8.30]; % Begin with GT measurement -end - - -% Compute variants and corresponding groupoid from euler angles, and -% convert the groupoids to Euler angles -[V,~] = YardleyVariants(ksi); -[G,~] = GroupoidVariantsCubic(V); -Ax = vector3d(G(:,2),G(:,3),G(:,4)); -G_eul = orientation('axis',Ax,'angle',G(:,1),CS_LT); - -% Classify groupoid as a misorientation in MTEX terms -rot = orientation('Euler',[0,0,0],CS_LT); -G_misos(1) = inv(rot) * rot; -G_misos(2:length(G_eul)+1) = inv(rot) * G_eul; -G_misos = transpose(G_misos); - -% Set tolerance (higher for KS and NW since expirimental usuallly varys -% from expected in those cases) -if (all(ksi==ksiKS) || all(ksi==ksiNW)) - tol = 1.25*degree; -else - tol = 0.85*degree; -end - -% make a copy of ebsd we can goof with -tmpEbsd = ebsd(ebsd.phaseId == 2); -tmpEbsd.phaseMap = [1,2,0]; -tmpEbsd.CSList = {tmpEbsd.CSList{1} tmpEbsd.CSList{2} tmpEbsd.CSList{end}}; - -% use the MTEX calcgrains to segment out Martensite grains -[grains,tmpEbsd.grainId] = calcGrains(tmpEbsd); - -% Isolate the grain boundaries and adjust the phase to singular -gB = grains.boundary; -gB.phase= 2; - -% Find the austenite grain boundaries based on the above tolerance -gB_aus = []; -for i = 1:length(G_misos) - if i == 1 - gB_aus = gB(angle(gB.misorientation,G_misos(i)) < tol); + + ksiKS = [5.26,10.30,10.53]; + ksiNW = [0,9.74,9.74]; + CS_LT = ebsd.CSList{2}; + CS_HT = ebsd.CSList{1}; + + % Allow the use of an initial Orientation Relationship guess. Otherwise, + % default to GT (ask eric, not sure what this is + if exist('OR_guess','var') + assert(isvector(OR_guess) && length(OR_guess) ==3 && isfloat(OR_guess),... + "OR_guess must be a length 3 array of floats") + ksi = OR_guess; % Begin with provided guess else - gB_aus = vertcat(gB_aus,gB(angle(gB.misorientation,G_misos(i)) < tol)); + ksi = [2.16,8.06,8.30]; % Begin with GT measurement end -end -% Merge grains together based on the chosen austenite grain boundaries -[mergedGrains,parentID]=merge(grains,gB_aus,'calcMeanorientation'); -mergedGrains.CS = CS_LT; - -% Find the maximum `austenite' grain and extract the Ebsd coordinates -% related to the original dataset -[~,maxGrnid] = max(mergedGrains.grainSize); -ParGrnId = find(parentID==maxGrnid); -EbId = []; -for i = 1:length(ParGrnId) - EbId = vertcat(EbId,find(tmpEbsd.grainId==ParGrnId(i))); -end - - -% Truncated ebsd dataset from `grain reconstruction' -tmpEbsd = tmpEbsd(EbId); -tmpEbsd.phase = 2; - -% Compute misorientation distribution function based on KS-OR and -% corresponding groupoid -hw = 2*degree; -psi=deLaValleePoussinKernel('halfwidth',hw); -MartModf = calcODF(G_misos,'kernel',psi); - -% Establish transformation matrices -[T2R,~]=calc_T2R(ksi,CS_HT,CS_LT); -[R2T,~]=calc_R2T(ksi,CS_HT,CS_LT); - -% From our EBSD region, transform the martensite to austenite and find the -% modal austenite orientation to use as our guess -aus_odf=calcODF(symmetrise(tmpEbsd.orientations)*T2R,'kernel',psi); -[~,AusOr] = max(aus_odf); - -% Call function to compute adjacency array and corresponding -% misorientations for the neighbgoring points -[adjpts,mori,~] = adjpt_moris(tmpEbsd,1); - -% In-plane and out-of-plane weights -modf_vals=eval(MartModf,mori); -modf_vals(modf_vals<0)=0; - -% If the input OR is KS or NW, adjust the parameterization for these. If -% it's GT or something similar, alter as well and include for twins. -if (all(ksi==ksiKS) || all(ksi==ksiNW)) - IP = [10,15]; - OP = [2,1.5e-1]; - % TwnWts = zeros(length(tmpEbsd),1); -else - IP = [7,8]; - OP = [2,1.5e-1]; - % Twns = vector3d([1,1,1;-1,-1,1;-1,1,1;1,-1,1]'); - % Twn_eul = orientation('axis',Twns,'angle',60*degree,ebsd{2}.CS); - % AusOrRot = rotation('euler',AusOr.phi1,AusOr.Phi,AusOr.phi2); - % TwnAus = AusOrRot * Twn_eul; - % TwnOrs = symmetrise(TwnAus)*R2T; - % TwnODF = calcODF(TwnOrs,'kernel',psi); - % TwnWts = (eval(TwnODF,tmpEbsd.orientations)+OP(1)).*OP(2); -end - -% Compute ODF for parent and establish the second set of out-of-plane -% weights -From_OrGuess=symmetrise(AusOr)*R2T; -From_GuessODF=calcODF(From_OrGuess,'kernel',psi); -OP2wts=(eval(From_GuessODF,tmpEbsd.orientations)+OP(1)).*OP(2); - -% Establish constant as first set of out-of-plane weights -% OP1wts = abs(1./OP2wts); -% OP1wts(OP1wts==Inf)=1e2; -% OP1wts = OP1wts + TwnWts; - -IPwts = (modf_vals+IP(1)).*IP(2); - -% Add source node -DiGraph=digraph; - -% Establish constant as first set of out-of-plane weights -OP1wts = abs(4./OP2wts); -OP1wts(OP1wts==Inf)=1e2; - -% Call function to set up the initial graph -[DiGraph,endnode,sinknode] = graph_setup(DiGraph,adjpts,IPwts,OP1wts,2); - -% Add edges corresponding to the guess reconstruction orientation -DiGraph=rmedge(DiGraph,(1:length(tmpEbsd))+endnode,sinknode); -DiGraph=addedge(DiGraph,(1:length(tmpEbsd))+endnode,sinknode,OP2wts); - -% Perform graph cut -[~,~,~,ct]=maxflow(DiGraph,1,sinknode); - -% Truncate ct and finalize the cut -ct(end) = []; -Cut = ct - endnode; - -% Convert the IDs in the grain to the Ids of the original scan and return -Aus_Grain_Ids = EbId(Cut); -% EbAus = tmpEbsd(Cut); -grain_remainder_Ids = EbId; -grain_remainder_Ids(Cut) = []; + + + % Compute variants and corresponding groupoid from euler angles, and + % convert the groupoids to Euler angles + [V,~] = YardleyVariants(ksi); + [G,~] = GroupoidVariantsCubic(V); + Ax = vector3d(G(:,2),G(:,3),G(:,4)); + G_eul = orientation('axis',Ax,'angle',G(:,1),CS_LT); + + % Classify groupoid as a misorientation in MTEX terms + rot = orientation('Euler',[0,0,0],CS_LT); + G_misos(1) = inv(rot) * rot; + G_misos(2:length(G_eul)+1) = inv(rot) * G_eul; + G_misos = transpose(G_misos); + + % Set tolerance (higher for KS and NW since expirimental usuallly varys + % from expected in those cases) + if (all(ksi==ksiKS) || all(ksi==ksiNW)) + tol = 1.25*degree; + else + tol = 0.85*degree; + end + + % make a copy of ebsd we can goof with + tmpEbsd = ebsd(ebsd.phaseId == 2); + tmpEbsd.phaseMap = [1,2,0]; + tmpEbsd.CSList = {tmpEbsd.CSList{1} tmpEbsd.CSList{2} tmpEbsd.CSList{end}}; + + % use the MTEX calcgrains to segment out Martensite grains + [grains,tmpEbsd.grainId] = calcGrains(tmpEbsd,'unitCell'); + + % Isolate the grain boundaries and adjust the phase to singular + gB = grains.boundary; + gB.phase= 2; + + % Find the austenite grain boundaries based on the above tolerance + gB_aus = []; + for i = 1:length(G_misos) + if i == 1 + gB_aus = gB(angle(gB.misorientation,G_misos(i)) < tol); + else + gB_aus = vertcat(gB_aus,gB(angle(gB.misorientation,G_misos(i)) < tol)); + end + end + % Merge grains together based on the chosen austenite grain boundaries + [mergedGrains,parentID]=merge(grains,gB_aus,'calcMeanorientation'); + mergedGrains.CS = CS_LT; + + % Find the maximum `austenite' grain and extract the Ebsd coordinates + % related to the original dataset + [~,maxGrnid] = max(mergedGrains.grainSize); + ParGrnId = find(parentID==maxGrnid); + EbId = []; + for i = 1:length(ParGrnId) + EbId = vertcat(EbId,find(tmpEbsd.grainId==ParGrnId(i))); + end + + + % Truncated ebsd dataset from `grain reconstruction' + tmpEbsd = tmpEbsd(EbId); + tmpEbsd.phase = 2; + + % Compute misorientation distribution function based on KS-OR and + % corresponding groupoid + hw = 2*degree; + psi=deLaValleePoussinKernel('halfwidth',hw); + MartModf = calcODF(G_misos,'kernel',psi); + + % Establish transformation matrices + [T2R,~]=calc_T2R(ksi,CS_HT,CS_LT); + [R2T,~]=calc_R2T(ksi,CS_HT,CS_LT); + + % From our EBSD region, transform the martensite to austenite and find the + % modal austenite orientation to use as our guess + aus_odf=calcODF(symmetrise(tmpEbsd.orientations)*T2R,'kernel',psi); + [~,AusOr] = max(aus_odf); + + % Call function to compute adjacency array and corresponding + % misorientations for the neighbgoring points + [adjpts,mori,~] = adjpt_moris(tmpEbsd,1); + + % In-plane and out-of-plane weights + modf_vals=eval(MartModf,mori); + modf_vals(modf_vals<0)=0; + + % If the input OR is KS or NW, adjust the parameterization for these. If + % it's GT or something similar, alter as well and include for twins. + if (all(ksi==ksiKS) || all(ksi==ksiNW)) + IP = [10,15]; + OP = [2,1.5e-1]; + % TwnWts = zeros(length(tmpEbsd),1); + else + IP = [7,8]; + OP = [2,1.5e-1]; + % Twns = vector3d([1,1,1;-1,-1,1;-1,1,1;1,-1,1]'); + % Twn_eul = orientation('axis',Twns,'angle',60*degree,ebsd{2}.CS); + % AusOrRot = rotation('euler',AusOr.phi1,AusOr.Phi,AusOr.phi2); + % TwnAus = AusOrRot * Twn_eul; + % TwnOrs = symmetrise(TwnAus)*R2T; + % TwnODF = calcODF(TwnOrs,'kernel',psi); + % TwnWts = (eval(TwnODF,tmpEbsd.orientations)+OP(1)).*OP(2); + end + + % Compute ODF for parent and establish the second set of out-of-plane + % weights + From_OrGuess=symmetrise(AusOr)*R2T; + From_GuessODF=calcODF(From_OrGuess,'kernel',psi); + OP2wts=(eval(From_GuessODF,tmpEbsd.orientations)+OP(1)).*OP(2); + + % Establish constant as first set of out-of-plane weights + % OP1wts = abs(1./OP2wts); + % OP1wts(OP1wts==Inf)=1e2; + % OP1wts = OP1wts + TwnWts; + + IPwts = (modf_vals+IP(1)).*IP(2); + + % Add source node + DiGraph=digraph; + + % Establish constant as first set of out-of-plane weights + OP1wts = abs(4./OP2wts); + OP1wts(OP1wts==Inf)=1e2; + + % Call function to set up the initial graph + [DiGraph,endnode,sinknode] = graph_setup(DiGraph,adjpts,IPwts,OP1wts,2); + + % Add edges corresponding to the guess reconstruction orientation + DiGraph=rmedge(DiGraph,(1:length(tmpEbsd))+endnode,sinknode); + DiGraph=addedge(DiGraph,(1:length(tmpEbsd))+endnode,sinknode,OP2wts); + + % Perform graph cut + [~,~,~,ct]=maxflow(DiGraph,1,sinknode); + + % Truncate ct and finalize the cut + ct(end) = []; + Cut = ct - endnode; + + % Convert the IDs in the grain to the Ids of the original scan and return + % EbAus = tmpEbsd(Cut); + grain_remainder_Ids = tmpEbsd.id; + grain_remainder_Ids(Cut) = []; + Aus_Grain_Ids = tmpEbsd.id(Cut); end \ No newline at end of file From 4b4b1927ee54e22b0b5aba9c40bb72573edc177c Mon Sep 17 00:00:00 2001 From: Brust Date: Fri, 12 Aug 2022 14:43:23 -0400 Subject: [PATCH 4/4] Script Updates Updated autoOR estimation to include unindexed points and fixed indexing error in find_twin_gbs --- Examples/Basics/Basics.m | 90 ++++++++++++++++++++--------- Mart2Aust/Core/AutoOR_estimation.m | 2 +- Mart2Aust/auxillary/find_twin_gbs.m | 11 ++-- 3 files changed, 70 insertions(+), 33 deletions(-) diff --git a/Examples/Basics/Basics.m b/Examples/Basics/Basics.m index d55d4d3..dacfe11 100644 --- a/Examples/Basics/Basics.m +++ b/Examples/Basics/Basics.m @@ -3,29 +3,19 @@ % Shows the simplest method for reconstructing the Austenite phase % from an EBSD scan of Martensite % ======================================================================= % - -% Create "AF96.ang" file from provided sample data if not already made -if size(dir('Af96.ang'),1) == 0 - data = h5read("../../Resources/EBSD/EBSD.hdf5","/AF96_small/AF001"); - header = h5readatt("../../Resources/EBSD/EBSD.hdf5","/AF96_small",'header'); - f = fopen('AF96.txt','w'); - fprintf(f,header); - fclose(f); - writematrix(data','AF96.txt','delimiter',' ','WriteMode','append') - movefile AF96.txt AF96.ang -end +tic % load some metadata about folder locations and do a flight check -Mart2Aust_Parent_folder = "C:\Users\agerlt\workspace\Mart2Aust"; +Mart2Aust_Parent_folder = "C:\Users\ALEBRUS\Desktop\Software\Matlab\Mart2Aust-main"; meta.Data_folder = Mart2Aust_Parent_folder + filesep +'Resources'+... filesep + 'EBSD'+ filesep +'AF96_small'; -meta.MTEX_folder = Mart2Aust_Parent_folder + filesep + 'MTEX'; +meta.MTEX_folder = Mart2Aust_Parent_folder + filesep + 'Mart2Aust' + filesep + 'MTEX'; meta.Functions_folder = Mart2Aust_Parent_folder + filesep + 'Mart2Aust'; meta.current_folder = string(pwd); meta.MTEX_Version = "mtex-5.7.0"; addpath(genpath(meta.Functions_folder)); -check_AusRecon_loaded; +% check_AusRecon_loaded; try disp(['Compatable MTEX version ', check_MTEX_version(meta.MTEX_Version), ' detected']) catch @@ -34,6 +24,8 @@ addpath(meta.MTEX_folder + filesep + meta.MTEX_Version); startup_mtex end +setMTEXpref('xAxisDirection','east') +setMTEXpref('zAxisDirection','outOfPlane') clear Aus_Recon_Parent_folders split_loc f data header Mart2Aust_Parent_folder @@ -45,40 +37,82 @@ % Load default reconstruction options and file location options = load_options("default"); -fname = dir('AF96.ang'); -% Use these two choices to generate some other location information -name = split(string(fname.name),'.'); -name = name(end-1); -location = [fname.folder, filesep, fname.name]; +pname = 'C:\Users\ALEBRUS\Documents\Microstructure_Work\Exxon_2022\Data'; +fname = '55-HR_EBSD1.ang'; +tname = fullfile(pname,fname); % Load a Martensite/Austenite Orientation Relationship and MDF halfwidth % that was previously measured on this scan (remove this line to recalculate % the values. See some of the other examples for more details on defining % and calculating the orientation relationship) -options.OR_ksi = [2.9057 7.7265 8.2255]; -options.OR_noise = 0.0281; +%% Alter Default Options +% options.OR_ksi = [2.7 9.2 9.4]; +% % options.OR_ksi = [2.4284 8.6749 8.8738]; +% % options.OR_ksi = [3.89,8.84,9.17]; +% options.OR_noise = 0.029; +% options.RGC_in_plane_m = 6; +% options.RGC_in_plane_b = 12; +% options.RGC_post_pre_m = 1.75e-1; +% options.RGC_post_pre_b = 0.6; +% options.OR_sampling_size = 2000; -% load the ebsd -original_ebsd = EBSD.load(location, ... +options.OR_plot_PAG_Mart_Guess = 1; +options.OR_plot_ODF_of_PAG_OR_guess = 1; + +%% Load EBSD and Format for Reconstruction +% Load original ebsd +original_ebsd = EBSD.load(tname, ... 'convertEuler2SpatialReferenceFrame','setting 2'); + +% Truncate ebsd if flagged +truncate = 1; +if truncate + ebsd = truncate_ebsd(original_ebsd,10,125,7,107); +else + ebsd = original_ebsd; +end + % Not all scans list phases in the same order: this command fixes that -ebsd = prep_for_Recon(original_ebsd,options); -clear original_ebsd +ebsd = prep_for_Recon(ebsd,options); + +% clear original_ebsd % at this point, we have identical scans with the following phase IDs: % 0 : Unindexed % 1 : Untransformed Parent (High temperature, or HT) % 2 : Transformed Child (Low Temperature, or LR) % 3 : Reconstructed Parent (starts empty) (Reconstructed, or R) +%% Auto OR Computation + % Determine the Orientation Relationship, if not already given -CS_HT =ebsd.CSList{1}; -CS_LT =ebsd.CSList{2}; +CS_HT = ebsd.CSList{1}; +CS_LT = ebsd.CSList{2}; [ebsd.opt.OR,ebsd.opt.HW,~] = AutoOR_estimation(ebsd,options); + +%% Calculate Misorientation Distribution Function + % Determine the Misorientation Distribution function for the LT phase [ebsd.opt.LT_MDF,ebsd.opt.psi] = calc_LT_MDF(CS_HT, CS_LT, ebsd.opt.HW, ebsd.opt.OR); +%% Perform Reconstruction +tic + +[Recon_ebsd, Recon_Likelihood] = Call_Reconstruction_Austin(ebsd,options); + +[a,AusGrns,E2627,E112_13] = CalcGrainSize(Recon_ebsd,3,1); + +Recon_Time = toc/60; +%% +tic +options.RGC_in_plane_m = 5; +options.RGC_in_plane_b = 3; +options.RGC_post_pre_m = 1.75e-1; +options.RGC_post_pre_b = 2; + % Perform the actual Reconstruction -Recon_ebsd = Call_Reconstruction(ebsd,options); +[Recon_ebsdOrig, Recon_LikelihoodOrig] = Call_Reconstruction(ebsd,options); +Recon_TimeOrig = toc/60; +%% Segment Variants % Segment the variants using the low temp and high temp maps variant_int_map = variants(ebsd, Recon_ebsd,options); diff --git a/Mart2Aust/Core/AutoOR_estimation.m b/Mart2Aust/Core/AutoOR_estimation.m index a6ed830..23a80b6 100644 --- a/Mart2Aust/Core/AutoOR_estimation.m +++ b/Mart2Aust/Core/AutoOR_estimation.m @@ -171,7 +171,7 @@ count = count+1; % Optimization function outside of ML add-on [MAPpars,~,~,~]=fminsearch(optimfunc,x0,MAP_options); - keyboard + if (MAPpars(1) < MAPpars(2) && MAPpars(3)) || count > 2 init_guess = 1; % Enforce constraint if need be diff --git a/Mart2Aust/auxillary/find_twin_gbs.m b/Mart2Aust/auxillary/find_twin_gbs.m index 875f814..6c25f69 100644 --- a/Mart2Aust/auxillary/find_twin_gbs.m +++ b/Mart2Aust/auxillary/find_twin_gbs.m @@ -34,12 +34,14 @@ % Merge parents and twins [Aus_gb_merged,ParentId] = merge(grains,All_Twin_Boundaries); - + % Loop through ebsd and replace grainIds with merged ones - for ii = 1:length(grains.id) - ebsd.grainId(ismember(ebsd.grainId,ii)) = ParentId(ii); + EbIds = zeros(length(ebsd),1); + for ii = 1:length(ParentId) + EbIds(ismember(ebsd.grainId,ii)) = ParentId(ii); end - + ebsd.grainId = EbIds; + % Assign parent orientations to the merged grain structure. for jj = 1:length(Aus_gb_merged.id) MergedGrns = grains.id(ismember(ParentId,jj)); @@ -47,6 +49,7 @@ ParOr = grains(MergedGrns(MaxId)).meanOrientation; Aus_gb_merged(jj).meanOrientation = ParOr; end + % NOTE 2: When you merged the grains, a lot of the grain orientations % were assigned values of nan. This fixes that, but under the % assumption that the largest surface area for each identified parent