gtDisorientation.m

function [mis_angle, mis_axis, mis_info] = gtDisorientation(ori1, ori2, symm, varargin)
% GTDISORIENTATION  Easy function to calculate the disorientation angle and axis
%                   between two sets of grains
%
%     [mis_angle, mis_axis, info] = gtDisorientation(ori1, ori2, symm, [varargin])
%     ----------------------------------------------------------------------------
%     INPUT:
%       ori1       = <double>   Orientations for grains set 1
%       ori2       = <double>   Orientations for grains set 2
%       symm       = <struct>   Symmetry operators list from gtCrystGetSymmetryOperators
%
%     OPTIONAL INPUT (varargin as a list of pairs, see parse_pv_pairs.m):
%       input      = <string>   Input orientation ({'rod'}, 'orimat', or 'euler')
%                                 rod:    input array sized 3xN
%                                 orimat: input array sized 3x3xN
%                                 euler:  input array sized 3xN
%       mode       = <string>   Mode of rotation ({'passive'} or 'active')
%       sort       = <string>   Sort disorientation axis components after taking
%                               absolute values (if empty, do nothing at all)
%                               ({''}, 'ascend' or 'descend')
%       latticepar = <double>   Lattice parameters (1x6 for hexagonal crystals)
%       convention = <string>   Convention used for the unit cell {'Y'}/'X' in
%                               case of hexagonal unit cell
%
%     OUTPUT:
%       mis_angle  = <double>   Misorientation angles between the 2 grain sets
%       mis_axis   = <double>   Misorientation axes between the 2 grain sets
%       mis_info   = <struct>   Extra info
%
%     Verison 002 24-10-2013 by LNervo, laura.nervo@esrf.fr
%       Added argument 'convention'
%
%     Version 001 27-06-2013 by YGuilhem, yoann.guilhem@esrf.fr

% Set default parameters and parse input parameters from varargin
par.input      = 'rod';
par.mode       = 'passive';
par.sort       = '';
par.latticepar = [];
par.convention = 'X';

[par, rej_pars] = parse_pv_pairs(par, varargin);

Symm = {symm.g3};

switch par.input
    case 'rod'
        Ngrain1 = size(ori1, 2);
        Ngrain2 = size(ori2, 2);
        oriConversion = @(x) gtMathsRod2OriMat(x);
    case 'orimat'
        Ngrain1 = size(ori1, 3);
        Ngrain2 = size(ori2, 3);
        oriConversion = @(x) x;
    case 'euler'
        Ngrain1 = size(ori1, 2);
        Ngrain2 = size(ori2, 2);
        oriConversion = @(x) gtMathsEuler2OriMat(x);
    otherwise
        gtError('gtDisorientation:wrong_input_type', ...
            'Allowed orientation types are ''rod'', ''orimat'', or ''euler''!');
end

if (Ngrain1 < 1 || Ngrain1 ~= Ngrain2)
    gtError('gtDisorientation:wrong_input_size', ...
        'This function needs equal number (N>1) of grains per set!');
end

switch par.mode
    case 'passive'
        calcDisorientation = @(x,y,z) sfPassiveDisorientation(x, y, z);
    case 'active'
        calcDisorientation = @(x,y,z) sfActiveDisorientation(x.', y.', z);
    otherwise
        gtError('gtDisorientation:wrong_rotation_mode', ...
            'Allowed rotation modes are only ''passive'' or ''active''!');
end

% Allocate output
mis_angle = zeros(Ngrain1, 1);

% Compute orientation matrices for sets 1 and 2
g1 = oriConversion(ori1);
g2 = oriConversion(ori2);

% Compute active/passive disorientation angle/axis
if (nargout > 1)
    mis_axis = zeros(Ngrain1, 3);
    
    tmp_info(Ngrain1).g1       = [];
    tmp_info(Ngrain1).g2       = [];
    tmp_info(Ngrain1).index_i  = [];
    tmp_info(Ngrain1).index_j  = [];
    tmp_info(Ngrain1).index_ij = [];
    tmp_info(Ngrain1).symm_i   = [];
    tmp_info(Ngrain1).symm_j   = [];
    tmp_info(Ngrain1).netg     = [];
    
    for ii = 1:Ngrain1;
        gA = g1(:, :, ii);
        gB = g2(:, :, ii);
        [mis_angle(ii), mis_axis(ii, :), mis_info] = calcDisorientation(gA, gB, Symm);
        
        tmp_info(ii) = gtAddMatFile(tmp_info(ii), mis_info, true);

        if (~isempty(par.latticepar)) && par.latticepar(6) == 120 % it must be hexagonal
            mis_axis_hex_cart(ii, 1:4) = gtCart2Hex(mis_axis(ii,:), par.latticepar, par.convention);
            tmp_info(ii).mis_axis_hex_cart = mis_axis_hex_cart(ii, :);

            if strcmp(par.mode, 'passive')
                [mis_axis_hex(ii, :), mis_axis_hex_int(ii, :)] = gtCrystMiller2FourIndexes(mis_axis(ii,:),'direction');
                tmp_info(ii).mis_axis_hex_int = mis_axis_hex_int(ii, :);
                tmp_info(ii).mis_axis_hex = mis_axis_hex(ii, :);
            end
        end
    end
    if nargout > 2
        mis_info = tmp_info;
    end
else
    for ii = 1:Ngrain1;
        gA = g1(:, :, ii);
        gB = g2(:, :, ii);
        mis_angle(ii) = calcDisorientation(gA, gB, Symm);
    end
end

if (~isempty(par.sort))
    mis_axis = sort(abs(mis_axis), 2, par.sort);
end

end % end of function

%% Sub-functions

function [theta, r, info] = sfPassiveDisorientation(gA, gB, Os)

N = uint32(length(Os));

% Compute all possible misorientation matrices
% This part could/should be optimized
netg = zeros(3, 3, 2*N^2);
DgAB = gA / gB; % equiv. and faster than gB * inv(gA);
for jj = 1:N
    Oj_DgAB = Os{jj} * DgAB;
    for ii = 1:N
        tmp = Os{ii} / Oj_DgAB; % equiv. and faster than Oj_DgAB * inv(Os{ii});
        netg(:, :, ii + N*(jj-1))   = tmp;
        netg(:, :, ii + N*(jj-1+N)) = tmp.';
        %netg(:, :,   (ii-1)*N + jj) = Os{ii} / Oj_DgAB; % equiv. and faster than Oj_DgAB * inv(Os{ii});
        %netg(:, :, (N+ii-1)*N + jj) = netg(:, :,   (ii-1)*N + jj).';
    end
end

% Convert them into angle/axis rotations
[all_theta, all_r] = gtMathsOriMat2AngleAxis(netg);

% Restrict choice to the fundamental zone
FZ_indices = find(all_r(1, :) >= all_r(2, :) & all_r(2, :) >= all_r(3, :) & all_r(3, :) >= 0);

% Get minimal rotation angle and axis
[theta, ind] = min(all_theta(:, FZ_indices));

% Convert FZ_index to all_index
index_ij = FZ_indices(ind);
r = all_r(:, index_ij).';

if nargout > 2
    info.g1 = gA;
    info.g2 = gB;
    [info.index_i, jj] = ind2sub([N 2*N], index_ij);
    info.index_j = jj;
    if info.index_j > N
        info.index_j = info.index_j - N;
    end
    info.index_ij = index_ij;
    info.symm_i = Os{info.index_i};
    info.symm_j = Os{info.index_j};
    info.netg = netg(:, :, FZ_indices(ind));
end
end

function [theta, r, info] = sfActiveDisorientation(gA, gB, Os)

N = length(Os);

% Compute all possible misorientation matrices
netg = zeros(3, 3, N);
for ii = 1:length(Os)
    netg(:, :, ii) = gA / (gB * Os{ii}); % equiv. and faster than gB * O{ii} * inv(gA)
end

% Convert them into angle/axis rotations
[all_theta, all_r] = gtMathsOriMat2AngleAxis(netg);

% Get minimal rotation angle and axis
[theta, ind] = min(all_theta(:));
r = all_r(:, ind).';

if nargout > 2
    info.g1 = gA;
    info.g2 = gB;
    info.index_i  = ind;
    info.index_j  = [];
    info.index_ij = [];
    info.symm_i = Os{ind};
    info.symm_j = [];
    info.netg = netg(:, :, ind);
end
end