gtVectorCryst2Lab.m

function Vs = gtVectorCryst2Lab(Vc, g, mode)
% GTVECTORLAB2CRYST  Express crystal row vector(s) in laboratory CS.
%
%     Vs = gtVectorCryst2Lab(Vc, g[, mode])
%     -----------------------------------
%
%     INPUT:
%       Vc   = <double>  Vector(s) expressed in the crystal CS.
%       g    = <double>  Array containing crystal orientation matrix(ces).
%       mode = <string>  Can be either 'diag' or a 3 coord vector {[3 1 2]},
%              <int>     This affect the output shape when multiple input 
%                        vectors and matrices, i.e. M>1 and N>1
%                        See further explanations about the output.
%
%     OUTPUT:
%       Vs   = <double>  Input vector(s) Vc expressed in the laboratory CS.
%                        Shape of output depends of input size (see below).
%
%     Note:
%       This function can be used in different modes:
%       - Vc = 1x3 array (1 row vector)  , g = 3x3xN array
%         output: N vectors expressed in the laboratory CS
%         output size: [N 3]
%
%       - Vc = Mx3 array (M row vectors) , g = 3x3 array
%         output: M vectors expressed in the laboratory CS
%         output size: [M 3]
%
%       - Vc = Mx3 array (M row vectors) , g = 3x3xN array
%         output: MxN vectors expressed in the laboratory CS
%         output size: [M 3 N] by default, can be modified using numeric mode
%           1: input vector coordinate dimension  (3)
%           2: input orientation matrix dimension (N)
%           3: input vector index dimension       (M)
%
%       - diag mode (only if M=N, need 'diag' mode, otherwise previous case)
%         Vc = Mx3 array (M row vectors) , g = 3x3xN array (M=N)
%         output: M vectors expressed in the laboratory CS
%         output size: [M 3 M]
%
%     Version 001 02-04-2013 by YGuilhem

% Get number of vectors (M) and orientation matrices (N)
M = size(Vc, 1);
N = size(g, 3);

diag = false; 
if ~exist('mode', 'var') || isempty(mode)
    mode = [3 1 2];
elseif isnumeric(mode)
    mode = mode(:).';
    if size(mode, 2) ~= 3 || ~all(between(mode, ones(size(mode)), 3*ones(size(mode))))
        gtError('gtVectorCryst2Lab:wrong_mode', ...
            'Numeric mode should have 3 coordinates in [1 3]!');
    end
elseif ischar(mode) && strcmp(mode, 'diag')
    if M ~= N
        gtError('gtVectorCryst2Lab:wrong_mode', ...
            '''diag'' mode needs as many vectors as matrices!');
    end
    diag = true;
else
    gtError('gtVectorCryst2Lab:wrong_mode', ...
        ['Unknown mode ''' mode '''!']);
end

% Check input orientation matrices g size
if (~isreal(g))
    gtError('gtVectorCryst2Lab:wrong_input_type', ...
        'Wrong input g type, wait a 3x3xN real array!');
elseif (size(g, 1) ~= 3 || size(g, 2) ~= 3)
    gtError('gtVectorCryst2Lab:wrong_input_size', ...
        'Input array g should be sized 3x3xN!');
end

% Check input vectors size
if size(Vc, 2) ~= 3
    gtError('gtVectorCryst2Lab:wrong_input_size', ...
        'Input array Vc should be sized Mx3!');
end

if M == 1 && N > 1
    % Transpose and reshape Vc
    Vc = reshape(Vc.', [3 1 M]);

    % Multiply and sum along 1st axis to reproduce Vs = Vc . g
    %Vs = squeeze(sum(bsxfun(@times, Vc, g), 1));
    Vs = squeeze(sum(Vc(:, [1 1 1], ones(N, 1)).*g, 1)).';

elseif M >= 1 && N == 1
%    % Transpose g to get gT, reshape it and expand it along 1st axis
%    gT = reshape(g.', [1 3 3]);
%
%    % Reshape Vc
%    Vc = reshape(Vc, [M 1 3]);
%
%    % Multiply and sum along 3rd axis to reproduce Vs = gT . Vc
%    Vs = sum(bsxfun(@times, gT, Vc), 3);

    % Can be simply done with this formulae 
    Vs = Vc * g;

elseif M == N && diag
    % Transpose and reshape Vc
    Vc = reshape(Vc.', [3 1 M]);

    % Multiply and sum along 1st axis to reproduce Vs = Vc . g
    %Vs = squeeze(sum(bsxfun(@times, Vc, g), 1)).';
    Vs = squeeze(sum(Vc(:, [1 1 1], :).*g, 1)).';

elseif M > 1 && N > 1
    % Transpose and reshape Vc
    Vc = reshape(Vc.', [3 1 1 M]);

    % Multiply and sum along 1st axis to reproduce Vs = Vc . g
    Vs = permute(squeeze(sum(bsxfun(@times, Vc, g), 1)), mode);

else
    gtError('gtVectorCryst2Lab:wrong_input_size', ...
            ['Cannot handle the case where size(Vc)=' num2str(size(Vc)) ...
            ' and length(g)=' num2str(N)]);
end

end % end of function