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