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, ...
'accelerate', false, ...
'mode', 'active', ... % ['active'], 'passive' : the mode of gtDisorientation
'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])
if conf.accelerate
disp('accelerated computation:')
[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 = [single(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) + 1);
for tmp_i = 1:numel(misorient_sigma)
indxs_sigma{tmp_i} = misorient_sigma(tmp_i).indices;
end
indxs_sigma{numel(misorient_sigma) + 1} = size(annealing_twins_Quaternion_origin, 2);
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)')];
misorient_sigma = [misorient_sigma, [1; 0; 0; 0]];
else
[annealing_twins_Quaternion_origin, twin_variant_path.paths] = 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 = [single(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
end
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
nof_neighboring_grains = numel(Neighboring_Grain_IDs); % number of neighboring grains
if conf.accelerate
%%%%%% The accelerated computations to determine the twin variants
% Find the twins among the neighboring grains
% calculate the disorientations between the current
% grain and its neighboring grains.
if strcmpi(conf.mode, 'active')
[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, 'mode', conf.mode);
mis_ax = mis_ax * gtMathsQuat2OriMat(Quaternions(:, gr_id_iter))';
mis_ax = sort(abs(mis_ax), 2, 'descend');
else
[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, 'mode', conf.mode, 'sort', 'descend');
end
% 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', 'mode', conf.mode);
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);
if ~isempty(Twin_Grain_IDs)
% 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', 'mode', conf.mode);
[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
end
else
mis_quat = gtMathsQuatProduct(...
Quaternions(:, gr_id_iter), ...
Quaternions(:, Neighboring_Grain_IDs), true);
min_misorients = zeros(1, nof_neighboring_grains);
inds_min_misorient = zeros(1, nof_neighboring_grains);
for tmp_i = 1:nof_neighboring_grains
tmp_misorients = gtDisorientation(...
annealing_twins_Quaternion_origin, ...
mis_quat(:, tmp_i) * ones(1, size(annealing_twins_Quaternion_origin, 2)), symm, 'input', 'quat', 'mode', conf.mode);
[min_misorients(tmp_i), inds_min_misorient(tmp_i)] = min(abs(tmp_misorients));
end
min_misorients = (min_misorients < misorient_tol);
inds_min_misorient = inds_min_misorient(min_misorients);
Twin_Grain_IDs = Neighboring_Grain_IDs(min_misorients);
end
% Check if the current grain has twins.
if ~isempty(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
% 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