gtRotateVolume.m

function vol = gtRotateVolume(vol, rot_tensor_angle, rot_axis, rot_center, interpolation)
% function vol = gtRotateVolume(vol, rot_tensor_angle, rot_axis, rot_center)
    if (any([size(rot_tensor_angle, 1), size(rot_tensor_angle, 2)] == 1))
        rot_angle = rot_tensor_angle;
        if (~exist('rot_axis', 'var') || isempty(rot_axis))
            rot_axis = [0 0 1];
        end
        rot_tensor = eye(3);
        for ii = 1:numel(rot_angle)
            rotcomp = gtMathsRotationMatrixComp(rot_axis(ii, :)', 'col');
            rot_tensor = gtMathsRotationTensor(rot_angle(ii), rotcomp) * rot_tensor;
%             rot_tensor = rot_tensor * gtMathsRotationTensor(rot_angle(ii), rotcomp);
        end
    else
        rot_tensor = rot_tensor_angle;
    end
    [vol_size(1), vol_size(2), vol_size(3)] = size(vol);
    if (~exist('rot_center', 'var') || isempty(rot_center))
        rot_center = (vol_size + 1) / 2;
    end
    if (~exist('interpolation', 'var') || isempty('interpolation'))
        interpolation = 'linear';
    end
    T1 = eye(4);
    T1(4, 1:3) = -rot_center;

    T2 = eye(4);
    T2(1:3, 1:3) = rot_tensor;
%     T2(1:3, 1:3) = rot_tensor';

    T3 = eye(4);
    T3(4, 1:3) = +rot_center;

    T = T1 * T2 * T3;

%     tform = affine3d(T);
%     vol = imwarp(vol, tform, 'OutputView', imref3d(size(vol)));

    % Old version. NOTE that rot_tensor should be un-transposed if this
    % version was to be used.
    tform = maketform('affine', T);
    R = makeresampler(interpolation, 'fill');
    vol = tformarray(vol, tform, R, [1 2 3], [1 2 3], vol_size, [], 0);
end