function W = rotation(V,alpha,beta,gamma)

% Compute mean of vertices for vertex list V:
center = mean(V)

% Initialize W:
nVertices = size(V,1)
W = zeros(nVertices,3);

% Translation matrix (translates the model to the point (0,0,0)):
T = [eye(3,3), -center'; zeros(1,3), 1];

% Translation back to the center of V:
Tback = [eye(3,3), center'; zeros(1,3), 1];

% Rotation matrices in homegeneuous coordinates:
Rx = [1 0 0 0; 0 cos(alpha) -sin(alpha) 0; 0 sin(alpha) cos(alpha) 0; 0 0 0 1]
Ry = [cos(beta) 0 sin(beta) 0;  0 1 0 0;  -sin(beta)  0  cos(beta) 0; 0 0 0 1]
Rz = [cos(gamma) -sin(gamma) 0 0; sin(gamma) cos(gamma) 0 0; 0 0 1 0; 0 0 0 1]

% Overall transformation matrix:
G = Tback * Rz * Ry * Rx * T;

% Homogeneous coordinates of V:
Vh = [V,ones(nVertices,1)];

% Rotate and translate the vertices:
Wh_t = G * Vh';
Wh = Wh_t';

% Back transform from homogenous to 3D coordinates:
W = Wh(:,1:3);