Perform the 3 axis rotation on a set of points.

rotationMatrix.m

function V = rotationMatrix(P, x, y, z)
    % rotationMatrix performs the 3 axis rotation on a
    % matrix P, which is a 3-by-p set of points, where
    % d is equal to 3 and represents the dimensionality
    % of the data and p represents the total number of
    % observations.
    %
    % x represents the amount of rotation along the x axis in radians
    % y represents the amount of rotation along the y axis in radians
    % z represents the amount of rotation along the z axis in radians
    %
    [d n] = size(P);
    V = zeros([d n]);

    M = [cos(y)*cos(z) -cos(x)*sin(z)+sin(z)*sin(y)*cos(z)  sin(x)*sin(z)+cos(x)*sin(y)*cos(z) ;
         cos(y)*sin(z)  cos(x)*cos(z)+sin(z)*sin(y)*sin(z) -sin(x)*cos(z)+cos(x)*sin(y)*sin(z) ;
         -sin(y)        sin(x)*cos(y)                       cos(x)*cos(y)                      ];

    if d == 3
        for i = 1:n
            V(:, i) = M * P(:, i);
        end
    end
end