function [common_refl_pl_table, common_refl_pl_celllist] = gtTwinCalcCommonReflections(parent_gr, twin_variants, hklsp, varargin)
% GTTWINCALCCOMMONREFLECTIONS  Function to calculate the common reflections
%                   between parent grain and twin variants.
%
%     common_refls = gtTwinCalcCommonReflections(twin_variants, parent_gr, tol, input_type, phase_id, output_sparse_flag)
%     ----------------------------------------------------------------------------
%     INPUT:
%       twin_variants       = <double>   Orientations for twin variants.
%                                        Default values are sigma 3 and 9
%                                        twin variants of FCC structure
%       parent_gr           = <double>   The orientation of the parent
%                                        grain. Default is the origin of
%                                        orientation space.
%     OPTIONAL INPUT:
%       'tol'               = <double>   The tolerance to check whether
%                                        hkls are equal.
%                                        Default value is 0.01.
%
%       'input_type'        = <string>   Input orientation ({'rod'}, 'orimat', or 'quat')
%                                        rod:    input array sized 3xN
%                                        orimat: input array sized 3x3xN
%                                        quat:  input array sized 4xN (default)
%                                        gr_id:  grain ids 1xN
%       'phase_id'          = <double>   phase ID. Default value is 1.
%       'output_sparse'     = <logical>  Output uses sparse (true) or not
%                                        (false). Default value is false.
%
%     OUTPUT:
%       common_refls  =  A table indicates common reflections. Its rows (vertical axis)
%                        represent the plane families of parent grain and
%                        its columns (horizontal axis) represent twin variants.
%                     e.g. common_refls(1, 2) = 3 means parameter.cryst.hklsp(:, 1)
%                          of parent grain and parameter.cryst.hklsp(:, 3) of the
%                          2nd twin variant have common reflections.
%
%     Verison 001 23-02-2021 by ZLiu, zheheng.liu@esrf.fr

    %% set default parameters
    if nargin < 2 || isempty(twin_variants)
        twin_variants = gtTwinCalcSigma3nVariantsQuat(2); % sigma 3n twin variants of FCC structure in crystal coordinates
    end
    conf = struct( ...
        'tol', 0.01, ...
        'input_type', 'quat', ...
        'phase_id', 1, ...
        'output_sparse', false);
    conf = parse_pv_pairs(conf, varargin);
    switch lower(conf.input_type)
        case 'rod'
            oriConversion = @(x) gtMathsRod2OriMat(x);
        case 'orimat'
            oriConversion = @(x) x;
        case 'gr_id'
            sample = GtSample.loadFromFile;
            oriConversion = @(x) gtMathsRod2OriMat(sample.phases{conf.phase_id}.R_vector(x, :)');
        case 'quat'
            oriConversion = @(x) gtMathsQuat2OriMat(x);
        otherwise
            gtError('gtTwinCalcCommonReflections:wrong_input_type', ...
                'Allowed orientation types are ''rod'', ''orimat'', or ''quat''!');
    end

    %% Calculate rotation matrix from twin variants to parent grain
    if ~nargin || isempty(parent_gr)
        g_twin_variants = permute(oriConversion(twin_variants), [2, 1, 3]);
    else
        g_twin_variants = gtMathsMatrixProduct(oriConversion(parent_gr), permute(oriConversion(twin_variants), [2, 1, 3]));
    end
    if ~exist('hklsp', 'var') || isempty(hklsp)
        parameters = gtLoadParameters();
        hklsp = parameters.cryst.hklsp;
    elseif isstruct(hklsp)
        if isfield(hklsp, 'hklsp')
            hklsp = hklsp.hklsp;
        elseif isfield(hklsp, 'cryst')
            hklsp = hklsp.cryst.hklsp;
        else
            error('Wrong input of HKL information')
        end
    end
    Parent_refl_ori = hklsp;
    % Transform the hkl of twin variants into coordinates of the parent.
    Twin_refl_ori = gtMathsMatrixProduct(g_twin_variants, Parent_refl_ori);

    %% Check whether the orientations are parallel
    % Compare each hkl_i of parent grain with each hkl_kj = R_k*hkl_j of
    % twin grains.
    % (||n_i - m_j||^2 < tol or ||n_i + m_j||^2 < tol) is simplified as
    % ||n_i||^2 + ||n_j||^2 - 2*|n_i'*n_j| < tol
    tmp_innerprod = 2 * abs(gtMathsMatrixProduct(permute(Parent_refl_ori, [2, 1, 3]), Twin_refl_ori));
    tmp_sqsum_t = sum(Twin_refl_ori.^2, 1);
    tmp_sqsum_p = permute(sum(Parent_refl_ori.^2, 1), [2, 1, 3]);
    tmp_sum = bsxfun(@plus, tmp_sqsum_p, tmp_sqsum_t);
    Common_table = bsxfun(@minus, tmp_sum, tmp_innerprod) < conf.tol;

    %% Transform the 3D logical matrix into 2D matrix
    % In the 2nd dimension of the 3D logical matrix "Common_table", each
    % column contains maximally one nonzero element. So we can represent
    % these columns by the indices of nonzero elements.
    if conf.output_sparse % use full or sparse matrix
        [tmp_ind1, tmp_ind_tw, tmp_ind2] = ind2sub(size(Common_table), find(Common_table));
        common_refl_pl_table = sparse(tmp_ind1, tmp_ind2, tmp_ind_tw);
    else
        common_refl_pl_table = uint8(sum(bsxfun(@times, Common_table, 1:size(Parent_refl_ori, 2)), 2));
        common_refl_pl_table = reshape(common_refl_pl_table, size(Parent_refl_ori, 2), size(g_twin_variants, 3));
    end

    % Extra output
    if nargout > 1
        common_refl_pl_celllist = cell(1, size(common_refl_pl_table, 1));
        hklind_shared = find(any(common_refl_pl_table, 2));
        tmpstruct = struct('grain', 0, 'HKL', [0, 0, 0]);
        for ii_hkl = 1:numel(hklind_shared)
            grainind_shared = find(common_refl_pl_table(hklind_shared(ii_hkl), :));
            common_refl_pl_celllist{hklind_shared(ii_hkl)} = tmpstruct(ones(1, numel(grainind_shared)));
            common_refl_pl_celllist{hklind_shared(ii_hkl)}(1).HKL = hklsp(:, hklind_shared(ii_hkl))';
            for ii_g = 1:numel(grainind_shared)
                common_refl_pl_celllist{hklind_shared(ii_hkl)}(ii_g+1) = struct(...
                    'grain', grainind_shared(ii_g), ...
                    'HKL', hklsp(:, common_refl_pl_table(hklind_shared(ii_hkl), grainind_shared(ii_g)))');
            end
        end
    end
end