Code covered by the BSD License  

Highlights from
SpinConv

SpinConv

by

 

Conversion from a rotation representation type to another

testSpinConv(testtype, verbose)
function testSpinConv(testtype, verbose)
%TESTSPINCONV   Comparing SpinConv with SpinCalc.
%
%    TESTSPINCONV(TESTTYPE) is equivalent to TESTSPINCONV(TESTTYPE, FALSE)
%
%    TESTSPINCONV(TESTTYPE, VERBOSE) compares either the output of SpinConv
%    and SpinCalc (if TESTTYPE == 'output') or their execution speed (if
%    TESTTYPE == 'speed').
%
%    76 different kinds of rotations (about x, y, z, and their
%    compositions) are converted from any kind of supported rotation
%    representation (DCM, EV, Q, EA121, EA232, EA313, EA131, EA212, EA323,
%    EA123, EA231, EA312, EA132, EA213, EA321) to any other kind. 
%
%    All the equations used by SpinConv and SpinCalc are tested, including
%    the four alternative kinds of equations used to convert from DCM to
%    quaternion. Both arrays with N=1 (single rotation) and arrays with N>1
%    (multiple rotations) are tested.
%
%    If VERBOSE == FALSE (default value), negligible differences smaller
%    than 1000 * EPS are not shown. If VERBOSE == TRUE, all differences are
%    shown. 
%
%    If TESTTYPE == 'speed', VERBOSE is ignored, N is increased up to 10000
%    by replicating and tiling the above mentioned set of rotations, and
%    only a selected subset of the above mentioned comparisons is
%    performed, to avoid redundancy.
%
%    Note: 
%        Output differences larger than 1000 * EPS are due to the fact
%        that, when SpinCalc performs a conversion to EV (Axis-angle
%        representation), in some cases it may return an angle MU larger
%        than 180 degrees, while SpinConv in that case equivalently
%        represents the same rotation as a rotation by an angle smaller
%        than 180 degrees (MU = 360 - MU) about an Euler vector pointing in
%        the opposite direction.
%        Output differences smaller than 1000 * EPS are due to the
%        different methods used by SpinConv and SpinCalc to convert angles
%        from degrees to radians and vice-versa. SpinConv uses the same
%        method used by MATLAB functions DEG2RAD and RAD2DEG (see
%        README.TXT). 
%
%    See also SpinConv, SpinCalc.

% 2013 February 25
% Paolo de Leva
% University "Foro Italico" 
% Rome, Italy
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Narg = nargin;
switch Narg
    case 0, testtype = 'output'; verbose=false;
    case 1, verbose=false;
    case 2
    otherwise, error( nargchk(0, 2, Narg) );
end
testtype = lower(testtype);
switch testtype
    case {'output', 'speed'}
        % Do nothing
    otherwise
        error(['Invalid value for first input argument ' ...
               '(use either ''ourput'' or ''speed'')'])
end

% Rotations about x, y, and z
R1 = rotations( pi/6);   % (3x3x3) by   30 degrees
R2 = rotations( pi*2/3); % (3x3x3) by  120 degrees
R3 = rotations(-pi/3);   % (3x3x3) by  -60 degrees
R4 = rotations(-pi*5/6); % (3x3x3) by -150 degrees
R = cat(3, R1, R2, R3, R4); % 3x3x12

% Permutation (x x x x   y y y y   z z z z)
neworder = [1 4 7 10  2 5 8 11  3 6 9 12];
R = R(:,:,neworder);
    
% Compositions of rotations about x, y, z
comb3 = combN([1 2 3 4], 3);
x = comb3(:,1);
y = comb3(:,2)+4;
z = comb3(:,3)+8;
for i = 4^3 : -1 : 1
    Rcomp(:,:,i) = R(:,:,z(i)) * R(:,:,y(i)) * R(:,:,x(i)); % 3x3x64
end
    
% Testing all kinds or rotations with all kinds of supported representations
clc
if strcmp(testtype, 'output')
    types  = {'DCM' 'EV' 'Q' ...
              'EA121' 'EA232' 'EA313' 'EA131' 'EA212' 'EA323' ...
              'EA123' 'EA231' 'EA312' 'EA132' 'EA213' 'EA321'};
    types2 = {'DCM' 'EV' 'Q' ...
              'EA123' 'EA231' 'EA312' 'EA132' 'EA213' 'EA321'};
    disp '1) Testing rotations about x'
    disp ' ', test_them([types2 {'EA313' 'EA212'}], R(:,:,1:4), testtype, verbose);
    disp ' '
    disp '2) Testing rotations about y'
    disp ' ', test_them([types2 {'EA121' 'EA323'}], R(:,:,5:8), testtype, verbose);
    disp ' '
    disp '3) Testing rotations about z'
    disp ' ', test_them([types2 {'EA232' 'EA131'}], R(:,:,9:12), testtype, verbose);
    disp ' '
    disp '4) Testing compositions (Rx*Ry*Rz)'
    disp ' ', test_them(types, Rcomp, testtype, verbose); % Compositions (Rx*Ry*Rz)
else
    clear R1 R2 R3 R4 R comb3 x y z
    idx = 1 : 64;
    idx = idx(ones(1,157), :);
    Rcomp = Rcomp(:,:,idx(:)');
    types  = {'DCM' 'EV' 'Q' 'EA121' 'EA123'};
    disp '1) Testing single rotation (N = 1)'
    disp ' '
    disp '   Conversion type         A             B           A-B           A/B'   
    test_it(types, Rcomp(:,:,1), testtype, verbose);
    disp ' ', disp '   A = SpinConv; B = SpinCalc; ms = milliseconds.'
    disp ' '
    disp '2) Testing multiple rotations (N = 10)'
    disp ' '
    disp '   Conversion type         A             B           A-B           A/B'   
    test_it(types, Rcomp(:,:,1:10), testtype, verbose);
    disp ' ', disp '   A = SpinConv; B = SpinCalc; ms = milliseconds.'
    disp ' '
    disp '3) Testing multiple rotations (N = 100)'
    disp ' '
    disp '   Conversion type         A             B           A-B           A/B'   
    test_it(types, Rcomp(:,:,1:100), testtype, verbose);
    disp ' ', disp '   A = SpinConv; B = SpinCalc; ms = milliseconds.'
    disp ' '
    disp '4) Testing multiple rotations (N = 1000)'
    disp ' '
    disp '   Conversion type         A             B           A-B           A/B'   
    test_it(types, Rcomp(:,:,1:1000), testtype, verbose);
    disp ' ', disp '   A = SpinConv; B = SpinCalc; ms = milliseconds.'
    disp ' '
    disp '5) Testing multiple rotations (N = 10000)'
    disp ' '
    disp '   Conversion type         A             B           A-B           A/B'   
    test_it(types, Rcomp(:,:,1:10000), testtype, verbose);
    disp ' ', disp '   A = SpinConv; B = SpinCalc; ms = milliseconds.'
end
    

function test_them(types, R, testtype, verbose)
% Test conversions from any supported representation type to any other.
    
    % STEP 1 - Single rotation (N=1)
    N = size(R, 3);
    verbose0 = false;
    for i = 1 : N
        test_it(types, R(:,:,i), testtype, verbose0);
    end

    % STEP 2 - Multiple rotations (N>1)
    test_it(types, R, testtype, verbose);


function test_it(types, R, testtype, verbose)
% Test conversions from any supported representation type to any other.

    % STEP 1 - Computing with SpinCalc all kinds of representations
    n = length(types);  
    rotation{1} = R;
    for j = n : -1 : 2
        conv = ['DCM' 'to' types{j}];
        rotation{j} = SpinCalc(conv, R, 10*eps, 1);
    end

    % STEP 2 - Comparing SpinConv with SpinCalc for all kinds of coversions
    for i = 1 : n
        for j = 1 : n
            if j~=i
                conv = [types{i} 'to' types{j}];
                switch testtype
                    case 'output', compare_output(conv, rotation{i}, verbose)
                    case 'speed',  compare_speed (conv, rotation{i})
                end
            end
        end
    end


function compare_speed(conv, rotation)

    % Comparing execution speed
    f0 = @() SpinCalc(conv, rotation, 10*eps, 0);
    f1 = @() SpinConv(conv, rotation, 10*eps, 0);
    t0 = timeit(f0);
    t1 = timeit(f1);
    ratio = 100 * (t1 / t0);
    delta = t1 - t0;
    blank = ' '; blanks = blank( ones(1,14-length(conv)) );
    conv = ['   ' conv blanks];
    disp ([conv sprintf('%11.3f ms',     t1*1000) ...
                sprintf('%11.3f ms',     t0*1000) ...
                sprintf('%+11.3f ms', delta*1000) ...
                sprintf('%10i', round(ratio)) '%'])


function compare_output(conv, rotation, verbose)

    % Comparing output
    out1 = SpinCalc(conv, rotation, 10*eps, 1);
    out2 = SpinConv(conv, rotation);
    diff = abs(out2-out1);
    maxdiff = max( diff(:) );
    if maxdiff > 1000*eps
        disp (['   ' conv ...
            ' - Diff. between SpinCalc and SpinConv: ' ...
            num2str(maxdiff) ' (see help)'])
    elseif verbose && maxdiff > eps
        disp (['   ' conv ...
            ' - Diff. between SpinCalc and SpinConv < 1000*EPS (see help)'])
    elseif verbose
        disp (['   ' conv ...
            ' - No diff. between SpinCalc and SpinConv'])
    end


function R = rotations(theta)
% Rotations of coordinate system relative to vector space, about x, y and z
% ("alias" transformations)
    s = sin(-theta);
    c = cos(-theta);
    Rx = [ 1  0  0;
           0  c -s;
           0  s  c ];
    Ry = [ c  0  s;
           0  1  0;
          -s  0  c ];
    Rz = [ c -s  0;
           s  c  0;
           0  0  1 ];
    R = cat(3, Rx, Ry, Rz);


function M = combN(V, N)
% N-compositions with repetition of elements taken from vector V
% See ALLCOMB, available on MATLAB Central File Exchange
% (http://www.mathworks.com/matlabcentral/fileexchange)

% List of all possible combinations of N elements
[Y{N:-1:1}] = ndgrid(V);
M = cat(N+1, Y{:});
% Reshape into matrix
M = reshape(M, [], N);

Contact us