gtDisorientation.m

function [mis_angle, mis_axis, 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] = gtDisorientation(ori1, ori2, [varargin])
%     -------------------------------------------------------------------------
%     INPUT:
%       ori1       = <double> Orientations for grains set 1
%       ori2       = <double> Orientations for grains set 2
%       symm       = <struct> Symmetry operators list from gtCrystGetSymmetryOperators
%       latticepar = <double> Lattice parameters (1x6 for hexagonal crystals)
%
%     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')
%
%     OUTPUT:
%       mis_angle  = <double> Misorientation angles between the 2 grain sets
%       mis_axis   = <double> Misorientation axes between the 2 grain sets
%
%     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.latticepar = [];
par = 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);
if nargout > 1
    if (~isempty(par.latticepar))
        mis_axis = zeros(Ngrain1, 4);
    else
        mis_axis = zeros(Ngrain1, 3);
    end
end

% Compute active/passive disorientation angle/axis
if nargout > 1
    tmp_info(Ngrain1).g1 = [];
    for ii = 1:Ngrain1;
        gA = g1(:, :, ii);
        gB = g2(:, :, ii);
        [mis_angle(ii), mis_axis(ii, :), tmp_info(ii)] = calcDisorientation(gA, gB, Symm);
    end
    if (~isempty(par.latticepar))
        mis_axis = gtCart2Hex(mis_axis, latticepar);
    end
    if nargout > 2
        info = tmp_info;
    end
else
    for ii = 1:Ngrain1;
        gA = g1(:, :, ii);
        gB = g2(:, :, ii);
        mis_angle(ii) = calcDisorientation(gA, gB, Symm);
    end
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, :));

% 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.index_ij = index_ij;
    [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.g1 = gA;
    info.g2 = gB;
    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.index_i = ind;
    info.g1 = gA;
    info.g2 = gB;
    info.symm = Os{ind};
    info.netg = netg(:, :, ind);
end
end