-
Notifications
You must be signed in to change notification settings - Fork 7
/
plotMap_blockWidths.m
241 lines (213 loc) · 8.35 KB
/
plotMap_blockWidths.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
function plotMap_blockWidths(job,varargin)
% % THIS SCRIPT WAS CONTRIBUTED BY: Dr Tuomo Nyyssönen
% % This script calculates the representative value for martensite block
% % widths by projecting all boundary points to the vector perpendicular
% % to the trace of the {111}a plane as per the following reference:
% % [S.Morito, H.Yoshida, T.Maki,X.Huang, Effect of block size on the
% % strength of lath martensite in low carbon steels, Mater. Sci. Eng.: A,
% % Volumes 438–440, 2006, Pages 237-240,
% % https://doi.org/10.1016/j.msea.2005.12.048]
%
% Syntax
% plotMap_blockWidths(job,varargin)
%
% Input
% job - @parentGrainreconstructor
% pGrainId - parent grain Id using the argument 'parentGrainId'
%
% Option
% noScalebar - Remove scalebar from maps
% noFrame - Remove frame around maps
if ~isempty(varargin) && any(strcmpi(varargin,'parentGrainId'))
pGrainId = varargin{find(strcmpi('parentGrainId',varargin)==1)+1};
else
warning('Argument "parentGrainId" (in single quotes) not specified. Block widths will be calculated for the EBSD map.');
pGrainId = job.parentGrains.id;
end
for ii = 1:length(pGrainId)
%% Define the parent grain
pGrain = job.parentGrains(job.parentGrains.id == pGrainId(ii));
% pEBSD = job.ebsd(pGrain);
% pEBSD = pEBSD(job.csParent);
%% Define the child grain(s)
clusterGrains = job.grainsPrior(job.mergeId == pGrainId(ii));
% cEBSD = job.ebsdPrior(job.ebsdPrior.id2ind(pEBSD.id));
% cEBSD = cEBSD(job.csChild);
%% Calculate martensite block widths
if length(pGrainId)==1
cGrains = clusterGrains(job.csChild);
[dBlock,zz,new_A_vec] = calcBlockWidth(job,pGrain,pGrainId(ii),cGrains);
else
cGrains{ii} = clusterGrains(job.csChild);
if ~isempty(cGrains{ii})
[dBlockNew{ii},~,~] = calcBlockWidth(job,pGrain,pGrainId(ii),cGrains{ii});
else
dBlockNew{ii} = nan;
end
%
% if job.isTransformed(find(job.mergeId == pGrainId(ii)))
% [dBlockNew{ii},~,~] = calcBlockWidth(job,pGrain,pGrainId(ii),cGrains{ii});
% else
% dBlockNew{ii} = nan;
% end
end
end
%% Define the text output format as Latex
setLabels2Latex
%% Define the window settings for a set of docked figures
% % Ref: https://au.mathworks.com/matlabcentral/answers/157355-grouping-figures-separately-into-windows-and-tabs
warning off
desktop = com.mathworks.mde.desk.MLDesktop.getInstance;
% % Define a unique group name for the dock using the function name
% % and the system timestamp
dockGroupName = ['plotMap_blockWidths_',char(datetime('now','Format','yyyyMMdd_HHmmSS'))];
desktop.setGroupDocked(dockGroupName,0);
bakWarn = warning('off','MATLAB:HandleGraphics:ObsoletedProperty:JavaFrame');
%% Plot the grains along with their traces and normals
drawnow;
figH = gobjects(1);
figH = figure('WindowStyle','docked');
set(get(handle(figH),'javaframe'),'GroupName',dockGroupName);
drawnow;
if length(pGrainId)==1
[~,mP] = plot(cGrains,dBlock,varargin{:});
hold all
ha(1) = quiver(squeeze(cGrains),squeeze(cross(zz,zvector)),'color',[1 0 0]);
ha(2) = quiver(squeeze(cGrains),squeeze(zz),'color',[0 1 0]);
ha(3) = quiver(squeeze(cGrains),squeeze(new_A_vec),'color', [0 0 1]);
legend(ha,'$111_a {\parallel} 011_m$ trace','$111_a {\parallel} 011_m$ normal','Mean of projected points')
else
for ii = 1:length(cGrains)
if ~isnan(dBlockNew{ii})
[~,mP] = plot(cGrains{ii},dBlockNew{ii},varargin{:});
end
hold all
end
% % https://au.mathworks.com/matlabcentral/answers/388193-how-do-i-convert-a-cell-array-of-different-size-cells-to-a-matrix
[dBlock,tf] = padcat(dBlockNew{:}); % concatenate, pad rows with NaNs
dBlock(~tf) = 0; % replace NaNs by zeros
dBlock = dBlock(:);
end
hold off
% Define the maximum number of color levels and plot the colorbar
colormap(flipud(bone));
caxis([0 round(max(dBlock))]);
colorbar('location','eastOutSide','lineWidth',1.25,'tickLength', 0.01,...
'YTick', [0:1:round(max(dBlock))],...
'YTickLabel',num2str([0:1:round(max(dBlock))]'), 'YLim', [0 round(max(dBlock))],...
'TickLabelInterpreter','latex','FontName','Helvetica','FontSize',14,'FontWeight','bold');
set(figH,'Name','Map: Trace & normal for block width calculation','NumberTitle','on');
if check_option(varargin,'noScalebar'), mP.micronBar.visible = 'off'; end
if check_option(varargin,'noFrame')
mP.ax.Box = 'off'; mP.ax.YAxis.Visible = 'off'; mP.ax.XAxis.Visible = 'off';
end
drawnow;
%% Plot the martensite block width histogram
drawnow;
figH = gobjects(1);
figH = figure('WindowStyle','docked');
set(get(handle(figH),'javaframe'),'GroupName',dockGroupName);
drawnow;
h = histogram(dBlock(dBlock>0),'Normalization', 'probability','faceColor',[162 20 47]./255);
set(gca,'FontSize',14);
xlabel('\bf Martensite block width [$\bf \mu$m]','FontSize',14,'FontWeight','bold');
ylabel('\bf Relative frequency [$\bf f$(g)]','FontSize',14);
set(figH,'Name','Histogram: Martensite block width','NumberTitle','on');
screenPrint('Step',['Figure ',num2str(figH.Number),': martensite block width histogram']);
drawnow;
% % Output histogram data in a table
class_range = h.BinEdges(2:end) - ((h.BinEdges(2)-h.BinEdges(1))/2);
disp(table(class_range',h.Values','VariableNames',{'blockWidth','Freq'}))
% % The figure and histogram show that block widths are consistently
% % smaller when calculated this way
%% Place first tabbed figure on top and return
warning on
allfigh = findall(0,'type','figure');
if length(allfigh) > 1
figure(length(allfigh)-1);
else
figure(1);
end
warning(bakWarn);
pause(1); % Reduce rendering errors
return
end
function [d_block_new,zz,new_A_vec] = calcBlockWidth(job,pGrain,pGrainId,cGrains)
%% Plot martensite block widths
% if check_option(varargin,'grains')
% Get all 111 vector3ds for each grain:
hh = Miller({1,1,1},{1,-1,1},{-1,1,1},{1,1,-1},job.p2c.CS);
hh = hh(job.packetId(cGrains.id))';
zz = pGrain.meanOrientation.project2FundamentalRegion.*hh;
% drawnow;
% figH = gobjects(1);
% figH = figure('WindowStyle','docked');
% set(get(handle(figH),'javaframe'),'GroupName',dockGroupName);
% drawnow;
% [~,mP] = plot(cGrains)
% hold all
% quiver(cGrains,cross(zz,zvector),'linecolor','r');
% quiver(cGrains,zz,'linecolor','g');
% hold off
% if check_option(varargin,'noScalebar'), mP.micronBar.visible = 'off'; end
% if check_option(varargin,'noFrame')
% mP.ax.Box = 'off'; mP.ax.YAxis.Visible = 'off'; mP.ax.XAxis.Visible = 'off';
% end
% drawnow;
% % Using the function at the bottom of this script, all grain boundary
% % points of a grain are projected to a vector going through the center of
% % that grain:
[p,~,~] = projectPoints2Vector(cGrains,rotate(cross(zz,zvector),...
rotation.byAxisAngle(zvector,90*degree)));
% % A representative value for the average halfwidth could be the mean of
% % the absolute values:
new_A = cellfun(@abs,p,'UniformOutput',false);
new_A = cellfun(@mean,new_A);
% % Vector form for visual verification:
new_A_vec = rotate(cross(zz,zvector),rotation.byAxisAngle(zvector,90*degree));
new_A_vec = normalize(new_A_vec).*new_A';
new_A_vec.antipodal = 1;
% % Calculate the block widths again and plot to visually verify
% % what the width perpendicular to the {111} trace looks like:
% % Calculate the assumed block width:
d_block_new = 2*new_A'.*sin(zz.theta);
end
function [p,px,py] = projectPoints2Vector(grains,v)
% This function projects all grain boundary points to a vector
% going through the center of that grain
%
% Syntax
% [p,px,py] = projectPoints2Vector(grains,v)
%
% Input:
% grains - @grain2d
% v - @vector3d, image plane components will be used
%
% Output:
% p - cell containing projection lengths for all boundary points
% px - cell containing projection point x coordinates
% py - cell containing projection point y coordinates
%
% Example:
% mtexdata testgrains
% v = caliper(grains(8),'shortest')
% p = projectPoints2Vector(grains(8),v)
if length(grains) ~= length(v)
error('number of input must be identical')
end
V = grains.V;
poly = grains.poly;
ce = grains.centroid;
for ii = 1:length(grains)
% Get vertices
Vg = V(poly{ii},:);
% Center vertices
Vg = [Vg(:,1)-ce(ii,1) Vg(:,2)-ce(ii,2)];
a = v(ii).y/v(ii).x;
b = -1;
c = Vg(:,2) - Vg(:,1)*a;
p{ii} = (a*Vg(:,2) - b*Vg(:,1))./sqrt(a^2+b^2);
px{ii} = (-a*c)/(a^2+b^2);
py{ii} = (-b*c)/(a^2+b^2);
end
end