Newer
Older
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
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
classdef GtTwinRelatedDomainRecognition < handle
properties (Access = public)
TRDs;
neighboring_TRDs;
twin_table;
neighboring;
phase_id = 1;
end
properties (Access = private)
annealing_cluster_type = 3;
end
methods (Access = public)
function self = GtTwinRelatedDomainRecognition(phaseID, misorient_tol, grain_ids, cryst, varargin)
% GtTwinRelatedDomainRecognition recognizes annealing twins from
% neighboring grains and labels the annealing twin boundaries.
%
% TRD_recognition =
% GtTwinRelatedDomainRecognition(phaseID, misorient_tol,
% grain_ids, cryst, varargin)
% --------------------------------------------------------------------------------------------------------
% INPUT:
% phaseID = <double> phase of the grains
% misorient_tol = <double> tolerance to determine the twins
% grain_ids = <double> Select the grains to calculate
% cryst = <struct> Structure containing
% crystal information.
% OPTIONAL INPUT:
% 'input' = <str> ['gr_id']: grain_ids contains grain IDs.
% 'rod': grain_ids contains Rodrigues
% vectors of the grains; if the
% 'neighbors' is not given, all the
% grains are considered to neighbor each
% other.
%
% 'neighbors' = neighboring table.
% 'neighbor_strategy' = <str> The method to create
% neighboring table, if
% neighboring table is not given.
% ['dil_vol']: based on
% dilated volume;
% 'bbox': based on bounding
% box.
% 'order_sigma' = <int> The largest inteval of the
% twinning order.
% 'pad' = <double> distance tolarence for neighboring
%
% OUTPUT:
% twin_clusters = <structure> A structure array contains twin
% clusters.
% The attribute "grain_id" is
% the grain ID of start grain;
% "twin_ids" is an array containing the
% IDs of the twins
% "twin_paths" is an array containing the
% labels of the twins
% neighboring_twin_clusters
% = <cell> Neighboring twin variants.
% twin_cands = <sparse> A table showing the twinning
% relationships. The label 1, 2, 3 and 4
% represent four sigma3 twin variants.
% The right of the label links the column
% grain and the left links the row grain.
% e.g. twin_cands(45, 67) = 123 represents
% the path from grain 45 to 67 is 1->2->3
% The decimal indicates the similar
% orientation.
% neighbors = <sparse> A table showing the neighboring
% relationships.
%
% Reference
% [1] Bryan W Reed, Roger W Minich, Robert E Rudd, and Mukul Kumar. The structure of the cubic coincident site lattice rotation group. Acta Crystallographica Section A: Foundations of Crystallography,60(3):263???277, 2004.
% [2] Cayron C. Multiple twinning in cubic crystals: geometric/algebraic study and its application for the identification of the ??3n grain boundaries[J]. Acta Crystallographica Section A: Foundations of Crystallography, 2007, 63(1): 11-29.
% [3] Heinz, A., & Neumann, P. (1991). Representation of orientation and disorientation data for cubic, hexagonal, tetragonal and orthorhombic crystals. Acta crystallographica section a: foundations of crystallography, 47(6), 780-789.
%% Input variables
if (~nargin || isempty(phaseID))
phaseID = 1;
end
self.phase_id = phaseID;
if (nargin < 2 || isempty(misorient_tol))
misorient_tol = 1.5; % tolerance
end
conf = struct('input', 'gr_id', ...
'neighbors', [], ...
'neighbor_strategy', 'dil_vol', ...
'order_sigma', 2, ...
'pad', 1);
conf = parse_pv_pairs(conf, varargin);
switch conf.input
case 'gr_id'
self.annealing_cluster_type = 3;
sample = GtSample.loadFromFile();
tot_grains = sample.phases{phaseID}.getNumberOfGrains; % the total number of grains
if (~exist('grain_ids', 'var') || isempty(grain_ids))
grain_ids = 1:tot_grains;
grain_ids = grain_ids(sample.phases{phaseID}.selectedGrains);
else
invalid_grain_ids = grain_ids > tot_grains | grain_ids < 1;
if (any(invalid_grain_ids))
error([mfilename ':wrong_argument'], 'grain_ids list contains %d invalid grain IDs.', sum(invalid_grain_ids))
disp('Invalid grain IDs:')
disp(grain_ids(invalid_grain_ids))
end
end
%% Load the quaternions of grains
Quaternions = gtMathsRod2Quat(sample.phases{phaseID}.R_vector.');
case 'rod'
self.annealing_cluster_type = [];
Quaternions = gtMathsRod2Quat(grain_ids);
tot_grains = size(grain_ids, 2);
grain_ids = 1:tot_grains;
conf.save = false;
conf.neighbors = true(tot_grains, tot_grains);
end
%% Find neighboring grains
if (size(conf.neighbors, 1) == tot_grains && size(conf.neighbors, 2) == tot_grains)
neighbors = logical(conf.neighbors(1:tot_grains, 1:tot_grains));
elseif ~isempty(conf.neighbors) && all(any(bsxfun(@eq, conf.neighbors(:), grain_ids(:)'), 2), 1)
neighbors = logical(sparse(tot_grains, tot_grains));
neighbors(conf.neighbors(:), conf.neighbors(:)) = true;
else
neighbors = gtFindNeighboringGrains(phaseID, grain_ids, conf.pad, sample, 'strategy', conf.neighbor_strategy); % the table indicates the neighboring grrains
end
%% Calculate the quaternions of all the possible annealing twin variants in crystal coordinates
if (~exist('cryst', 'var') || isempty(cryst))
param = gtLoadParameters();
cryst = param.cryst;
end
space_group = cryst(phaseID).spacegroup;
latticepar = cryst(phaseID).latticepar;
symm = cryst(phaseID).symm;
order_sigma = conf.order_sigma;
twin_variant_path = struct(...
'paths', [], 'paths_inv', [], ...
'base', 10 .^ (0:order_sigma-1));
fcc_space_group = {196, 202, 203, 209, 210, 216, 219, 225, 226, 227, 228};
bct_space_group = 123;
switch space_group
case fcc_space_group
cryst_type = 1;
case bct_space_group
cryst_type = 2;
otherwise
cryst_type = -1;
end
if any(cryst_type == [1, 2])
[annealing_twins_Quaternion_origin, twin_variant_path.paths, misorient_sigma] = gtTwinCalcSigma3nVariantsQuat(order_sigma, space_group, symm, latticepar); % sigma 3n twin variants of FCC structure in crystal coordinates
% include its own orientation with label 0
annealing_twins_Quaternion_origin = [annealing_twins_Quaternion_origin, [1; 0; 0; 0]];
twin_variant_path.paths = [twin_variant_path.paths, [1e-5; zeros(order_sigma-1, 1)]]; % 1e-5 denotes the similar orientation
twin_variant_path.paths_inv = twin_variant_path.paths;
% reverse the labels
for tmp_i = 1:size(twin_variant_path.paths, 2)
twin_variant_path.paths_inv(twin_variant_path.paths(:, tmp_i) > 0, tmp_i) = flipud(twin_variant_path.paths(twin_variant_path.paths(:, tmp_i) > 0, tmp_i));
end
indxs_sigma = cell(1, numel(misorient_sigma));
for tmp_i = 1:numel(misorient_sigma)
indxs_sigma{tmp_i} = misorient_sigma(tmp_i).indices;
end
misorient_sigma = [cosd(cat(2, misorient_sigma.rot_ang) / 2); ...
bsxfun(@times, sind(cat(2, misorient_sigma.rot_ang) / 2), cat(1, misorient_sigma.rot_axis)')];
else
error([mfilename ':unknown crystal structure'], 'the function does not currently support this crystal structure.');
return;
end
%% Find the twin variants
neighboring_twin_clusters = cell(0);
flag_TRD_graph_node = uint8(zeros(tot_grains, 1)); % 2: The central grain in a twin cluster; 1: a twin; 0: not twin
twin_cands = sparse(tot_grains, tot_grains); % For each grain, store an array containing the misorientations between its annealing twin and its neighboring grains.
% some variables to simplify the computations.
tmp_dup_indxs = struct('twins', [], 'neighbors', []); % indices to duplicate and extend the quaternions to avoid redundant loop.
flags_unused_grs = true(tot_grains, 1);
indx_clusters = 0; % the index for twin clusters
queue_cluster = zeros(2 + tot_grains, 1);
% loop for the grains
gauge = GtGauge(numel(grain_ids), 'Calculating misorientations: ');
for tmp_grain_id_iter = 1:numel(grain_ids)
% If the grain has been computed, we skip this grain.
if flags_unused_grs(grain_ids(tmp_grain_id_iter))
% Initialize the queue of cluster members and push the grain
% into the queue as the first element.
queue_cluster(1:3) = [2; 3; grain_ids(tmp_grain_id_iter)];
% loop until queue is empty.
while(queue_cluster(1) ~= queue_cluster(2)) %(queue_cluster.queue_head ~= queue_cluster.queue_end)
% Pop the grain ID from the queue
queue_cluster(1) = queue_cluster(1) + 1;
gr_id_iter = queue_cluster(queue_cluster(1));
% As we are going to computing this grain, change its flag.
flags_unused_grs(gr_id_iter) = false;
% find the neighboring grains that have not been computed
Neighboring_Grain_IDs = and(neighbors(:, gr_id_iter), flags_unused_grs);
Neighboring_Grain_IDs = find(Neighboring_Grain_IDs);
if ~isempty(Neighboring_Grain_IDs) % make sure there exist neighboring grains that have not been computed.
% Check if the current grain is a twin of the previous grains. If so, add it
% into the current cluster.
if flag_TRD_graph_node(gr_id_iter)
indx_twins_current_cluster = indx_twins_current_cluster + 1;
neighboring_twin_clusters{indx_clusters}(indx_twins_current_cluster).grain_id = gr_id_iter;
end
%%%%%% The computations to determine the twin variants
%%%%%% are in the part I and II below.
%%%%%% Part I: Find the twins among the neighboring grains
nof_neighboring_grains = numel(Neighboring_Grain_IDs); % number of neighboring grains
% calculate the disorientations between the current
% grain and its neighboring grains.
[mis_ang, mis_ax] = gtDisorientation(...
Quaternions(:, gr_id_iter) * ones(1, nof_neighboring_grains), ...
Quaternions(:, Neighboring_Grain_IDs), symm, 'input', 'quat', 'constrain_to_sst', false, 'sort', 'descend');
% transform the rotation angles and axes of disorientations into quaternions
mis_quat = [cosd(mis_ang(:)/2), bsxfun(@times, sind(mis_ang(:)/2), mis_ax)]';
[tmp_dup_indxs.twins, tmp_dup_indxs.neighbors] = self.sfDuplicateIndices(size(misorient_sigma, 2), nof_neighboring_grains);
% determine what types of twin variants (sigma3, sigma9, sigma27a, sigma27b or none)
% the neighboring grains are in disorientation space[3].
tmp_misorients = gtDisorientation(...
misorient_sigma(:, tmp_dup_indxs.twins), ...
mis_quat(:, tmp_dup_indxs.neighbors), symm, 'input', 'quat');
tmp_misorients = reshape(tmp_misorients(:), [nof_neighboring_grains, size(misorient_sigma, 2)]);
[min_misorients, inds_min_misorient] = min(abs(tmp_misorients), [], 2);
% select the twins among the neighboring grains
min_misorients = (min_misorients < misorient_tol);
inds_min_misorient = inds_min_misorient(min_misorients);
Twin_Grain_IDs = Neighboring_Grain_IDs(min_misorients);
%%%%%% Part I end %%%%%%%
% Check if the current grain has twins.
if numel(Twin_Grain_IDs)
% Check if the current grain is a twin of any other grain. If not, set
% it as the central grain of the new cluster.
% Namely, if it has not been put into any twin cluster and
% it has twins, then it must belong to a new twin cluster.
% And we regard it as the central grain of this new cluster.
if ~flag_TRD_graph_node(gr_id_iter)
flag_TRD_graph_node(gr_id_iter) = 2; % root of the graph indicated by 2
indx_clusters = indx_clusters + 1; % new cluster
indx_twins_current_cluster = 1; % initialize the iteration index for twins
neighboring_twin_clusters{indx_clusters} = self.sfCreateClusterStruct(gr_id_iter);
end
%%%%%% Part II: Determine the specific twin variants
% Calculate the substraction between the rotations of
% the current grain and its twins.
mis_quat = gtMathsQuatProduct(...
Quaternions(:, gr_id_iter), ...
Quaternions(:, Twin_Grain_IDs), true);
% loop for its twins
for tmp_i = numel(Twin_Grain_IDs):-1:1
% For example, if the twin is sigma27a, then "tmp_twin_inds"
% is the indices of all the sigma27a twin variants.
tmp_twin_inds = indxs_sigma{inds_min_misorient(tmp_i)};
% Determine which specific twin variant the twin is among the twin variants at indices "tmp_twin_inds".
tmp_misorients = gtDisorientation(...
annealing_twins_Quaternion_origin(:, tmp_twin_inds), ...
mis_quat(:, tmp_i) * ones(1, numel(tmp_twin_inds)), symm, 'input', 'quat');
[min_misorients, inds_min] = min(abs(tmp_misorients));
if abs(min_misorients) < misorient_tol
% Get the index of the twin variant.
inds_min_misorient(tmp_i) = tmp_twin_inds(inds_min);
else
inds_min_misorient(tmp_i) = [];
mis_quat(:, tmp_i) = [];
Twin_Grain_IDs(tmp_i) = [];
end
end
%%%%%% Part II end %%%%%%
% put the labels in the table
twin_cands(Twin_Grain_IDs, gr_id_iter) = sum(bsxfun(@times, ...
double(twin_variant_path.paths(:, inds_min_misorient)'), ...
twin_variant_path.base), 2);
twin_cands(gr_id_iter, Twin_Grain_IDs) = sum(bsxfun(@times, ...
double(twin_variant_path.paths_inv(:, inds_min_misorient)), ...
twin_variant_path.base'), 1);
% put the twin IDs and labels in the cell
neighboring_twin_clusters{indx_clusters}(indx_twins_current_cluster).twin_ids = find(twin_cands(:, gr_id_iter) > 0).';
neighboring_twin_clusters{indx_clusters}(indx_twins_current_cluster).twin_paths = round(full(twin_cands(twin_cands(:, gr_id_iter) > 0, gr_id_iter)))';
% Update the quaternions of the twins to match the
% twin related domain. However, the updated value
% is just approximately equivalent to old value, so
% this step should be improved.
tmp_inds = ~flag_TRD_graph_node(Twin_Grain_IDs);
inds_min_misorient = inds_min_misorient(tmp_inds);
Twin_Grain_IDs = Twin_Grain_IDs(tmp_inds);
Quaternions(:, Twin_Grain_IDs) = gtMathsQuatProduct(...
Quaternions(:, gr_id_iter), ...
annealing_twins_Quaternion_origin(:, inds_min_misorient));
flag_TRD_graph_node(Twin_Grain_IDs) = 1; % branch of the graph indicated by 1
% extend the member queue of the cluster
queue_cluster((queue_cluster(2) + 1):(queue_cluster(2) + numel(Twin_Grain_IDs))) = Twin_Grain_IDs.';
queue_cluster(2) = queue_cluster(2) + numel(Twin_Grain_IDs);
end
end
gauge.incrementAndDisplay()
end
end
end
%% To sort out the twin clusters and store them in "twin_clusters"
twin_clusters = self.sfCreateClusterStruct([]);
gauge = GtGauge(length(neighboring_twin_clusters), 'Sorting out the twin clusters: ');
% loop of the clusters
for iter_cluster = 1:length(neighboring_twin_clusters)
gauge.incrementAndDisplay()
missing_twin_ids_cluster_bool = false(tot_grains, 1);
for tmp_i = length(neighboring_twin_clusters{iter_cluster}):-1:2
if isempty(neighboring_twin_clusters{iter_cluster}(tmp_i).twin_ids)
neighboring_twin_clusters{iter_cluster}(tmp_i) = [];
else
missing_twin_ids_cluster_bool(neighboring_twin_clusters{iter_cluster}(tmp_i).grain_id) = true;
missing_twin_ids_cluster_bool(neighboring_twin_clusters{iter_cluster}(tmp_i).twin_ids) = true;
end
end
twin_clusters(iter_cluster) = neighboring_twin_clusters{iter_cluster}(1);
missing_twin_ids_cluster_bool(twin_clusters(iter_cluster).grain_id) = false;
missing_twin_ids_cluster_bool(twin_clusters(iter_cluster).twin_ids) = false;
nof_existing_twins = numel(twin_clusters(iter_cluster).twin_ids);
nof_missing_twins = sum(missing_twin_ids_cluster_bool);
if (nof_missing_twins > 0) % check whether there is any missing twin
% assign sufficient space to store the missing twins
twin_clusters(iter_cluster).twin_ids = [twin_clusters(iter_cluster).twin_ids, zeros(1, nof_missing_twins)];
twin_clusters(iter_cluster).twin_paths = [twin_clusters(iter_cluster).twin_paths, zeros(1, nof_missing_twins)];
for twin_iter = 2:length(neighboring_twin_clusters{iter_cluster})
twin_ids_sub_cluster = neighboring_twin_clusters{iter_cluster}(twin_iter).twin_ids;
missing_twin_ids_sub_cluster_ind = find(missing_twin_ids_cluster_bool(twin_ids_sub_cluster));
missing_twin_ids_sub_cluster = twin_ids_sub_cluster(missing_twin_ids_sub_cluster_ind);
missing_twin_ids_cluster_bool(missing_twin_ids_sub_cluster) = false;
shared_twin_path = twin_clusters(iter_cluster).twin_paths(twin_clusters(iter_cluster).twin_ids == neighboring_twin_clusters{iter_cluster}(twin_iter).grain_id);
for tmp_i = 1:numel(missing_twin_ids_sub_cluster_ind)
% add the grain into the twin list
nof_existing_twins = nof_existing_twins + 1;
twin_clusters(iter_cluster).twin_ids(nof_existing_twins) = missing_twin_ids_sub_cluster(tmp_i);
% Calculate the new label
twin_clusters(iter_cluster).twin_paths(nof_existing_twins) = self.sfCombineLabels(shared_twin_path, ...
neighboring_twin_clusters{iter_cluster}(twin_iter).twin_paths(missing_twin_ids_sub_cluster_ind(tmp_i)), cryst_type);
end
end
end
end
% outputs
self.TRDs = twin_clusters;
self.neighboring_TRDs = neighboring_twin_clusters;
self.twin_table = twin_cands;
self.neighboring = neighbors;
end
function saveToSample(self, rspace_oversize)
if (~exist('rspace_oversize', 'var') || isempty(rspace_oversize))
rspace_oversize = 1;
end
phaseID = self.phase_id;
twin_clusters = self.TRDs;
% define the cluster_type of annealing twin clusters
annealing_cluster_flag = self.annealing_cluster_type;
sample = GtSample.loadFromFile;
sample.phases{phaseID} = GtPhase.loadobj(sample.phases{phaseID});
% the number of all the existing clusters in the sample file
nof_previous_clusters = length(sample.phases{phaseID}.clusters);
% add the annealing twin clusters into the class
for tmp_i = 1:length(twin_clusters)
% First row includes the grain IDs. The first grain is the
% central grain.
% Second row includes the corresponding labels.
incl_ids = [twin_clusters(tmp_i).grain_id, twin_clusters(tmp_i).twin_ids];
% bounding box in spatial space
% using the bounding boxes in sample.mat or loading grains to
% use gt6DMergeRealSpaceVolumes? Latter func involves det_index
bb_rs = max(sample.phases{phaseID}.boundingBox(:, 1:3) + sample.phases{phaseID}.boundingBox(:, 4:6) - 1, [], 1);
bb_rs = [min(sample.phases{phaseID}.boundingBox(:, 1:3), [], 1), bb_rs(:, 1:3)];
bb_rs = [bb_rs(:, 1:3), bb_rs(:, 4:6) - bb_rs(:, 1:3) + 1];
bb_rs_over_sz = ceil(((rspace_oversize - 1)/2) * bb_rs(:, 4:6));
bb_rs = bb_rs + [-bb_rs_over_sz, 2*bb_rs_over_sz];
% bounding boxes in orientation space
bb_o_space = [];
% GtPhase.makeCluster is to construct the cluster with given
% parameters
tmp_cluster = GtPhase.makeCluster(...
incl_ids, ...
sample.phases{phaseID}.getR_vector(incl_ids(1)), ...
bb_rs, ...
bb_o_space, ...
annealing_cluster_flag, ...
[0, twin_clusters(tmp_i).twin_paths], ...
false);
sample.phases{phaseID}.clusters(end + 1) = tmp_cluster;
end
% remove previous annealing twin clusters
for tmp_i = nof_previous_clusters:-1:1
if sample.phases{phaseID}.clusters(tmp_i).cluster_type == annealing_cluster_flag
sample.phases{phaseID}.clusters(tmp_i) = [];
end
end
% save into the file
sample.saveToFile;
end
end
methods (Access = protected)
% duplicate the indices
function [dup_inds_x, dup_inds_y] = sfDuplicateIndices(~, nof_x, nof_y)
% equivalent to [dup_inds_x, dup_inds_y] = meshgrid(1:nof_x, 1:nof_y).
dup_inds_x = ones(nof_y, 1) * (1:nof_x);
dup_inds_y = (1:nof_y)' * ones(1, nof_x);
dup_inds_x = dup_inds_x(:);
dup_inds_y = dup_inds_y(:);
end
% The structure representing the twin cluster.
function cluster_struct = sfCreateClusterStruct(~, gr_id)
cluster_struct = struct(...
'grain_id', gr_id, ...
'twin_ids', [], ...
'twin_paths', []);
end
% sub-function to combine the labels
function new_label = sfCombineLabels(~, label_first_half, label_second_half, cryst_type)
if label_first_half && label_second_half
label_first_half = num2str(label_first_half);
label_second_half = num2str(label_second_half);
% Remove all the pairs of repetitive labels.
if cryst_type == 1
opposite_func = @(x, y) x==y;
else
opposite_func = @(x, y) ((abs(x - y) == 1) & (abs(x + y - 5) == 2));
end
while(~isempty(label_second_half) && ~isempty(label_first_half) && opposite_func(str2double(label_second_half(end)), str2double(label_first_half(1))))
label_second_half(end) = [];
label_first_half(1) = [];
end
if isempty([label_second_half(:); label_first_half(:)])
new_label = 0; % represent the grain has the same orientation with the central grain
else
new_label = round(str2double([label_second_half(:); label_first_half(:)]));
end
else
new_label = label_first_half + label_second_half;
end
end
end
end