find common reflection HKL planes according to the Rodrigues vectors.

Signed-off-by: Zheheng Liu <zheheng.liu@esrf.fr>

zUtil_Twins/gtTwinCalcCommonReflections.m

+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
\ No newline at end of file