MATLAB rotation problem

Hey all, I need rotated versions of a 3D shepp logan phantom and it's corresponding rotation matrix. Now here's the thing, I use a function called phantom3d to create the 3D SLP, the function allows euler angles to specify a rotation. So for example:

phi = 45;
theta = 45;
psi = 45;

%just a matrix of inputs to create the shepp logan phantom
e =[  1  .6900  .920  .810      0       0       0    开发者_如何学C  0+phi      0+theta      0+psi
    -.8  .6624  .874  .780      0  -.0184       0      0+phi      0+theta      0+psi
    -.2  .1100  .310  .220    .22       0       0    -18+phi      0+theta     10+psi
    -.2  .1600  .410  .280   -.22       0       0     18+phi      0+theta     10+psi
     .1  .2100  .250  .410      0     .35    -.15      0+phi      0+theta      0+psi
     .1  .0460  .046  .050      0      .1     .25      0+phi      0+theta      0+psi
     .1  .0460  .046  .050      0     -.1     .25      0+phi      0+theta      0+psi
     .1  .0460  .023  .050   -.08   -.605       0      0+phi      0+theta      0+psi
     .1  .0230  .023  .020      0   -.606       0      0+phi      0+theta      0+psi
     .1  .0230  .046  .020    .06   -.605       0      0+phi      0+theta      0+psi   ];

img = phantom3d(e, 50);

Now according to the literature you can calculate a rotation matrix using:

phi = ((phi + 180)/180).*pi;
theta = (theta/180).*pi;
psi = (psi/180).*pi;

cphi = cos(phi);
sphi = sin(phi);
ctheta = cos(theta);
stheta = sin(theta);
cpsi = cos(psi);
spsi = sin(psi);

% Euler rotation matrix
alpha = [cpsi*cphi-ctheta*sphi*spsi   cpsi*sphi+ctheta*cphi*spsi  spsi*stheta;
        -spsi*cphi-ctheta*sphi*cpsi  -spsi*sphi+ctheta*cphi*cpsi cpsi*stheta;
        stheta*sphi  

However, if I compare the image i create using the phantom3d with a function which applies the rotation matrix on a non-rotated image, they don't rotate in the same way. The code to view the rotated versionof this image is:

img = phantom3d(50);
szout = size(img);
Cf = eye(4);
Cf(1:3, 4) = -szout/2;
Co = Cf;
%previously created alpha
alpha(4,4) = 1;
%Cf & Co are used for translations
Rmatrix = inv(Cf) * alpha * Co;
[x, y, z]=ndgrid(single(1:szout(1)), single(1:szout(2)), single(1:szout(3)));
xyz = [x(:) y(:) z(:) ones(numel(x),1)]*Rmatrix(1:3,:)';
xyz = reshape(xyz,[szout 3]);
img2 = interpn(single(img), xyz(:,:,:,1),xyz(:,:,:,2),xyz(:,:,:,3), 'cubic', 0);

So I actually need to have img & img2 to be the same, but it's not. I found some case where I set psi, phi & theta to 45 and then add 180 to phi when creating img2 it gives the same result, so there is some relation to it but I can't seem to find it.

Anyone have any ideas, suggestions, help?

Thx a lot

Problem solved, apparently the rotation over the x axis is different in this function. Normally when u calculate the rotation matrix for euler angles they state it to be:

D = [cos(phi) sin(phi) 0 -sin(phi) cos(phi) 0 0 0 1];

C = [1 0 0 0 cos(theta) sin(theta) 0 -sin(theta) cos(theta)];

B = [cos(psi) sin(psi) 0 -sin(psi) cos(psi) 0 0 0 1];

R = B*C*D;

however in my case the C was different, namely a rotation over the old y axis:

C = [cos(theta) 0 -sin(theta) 0 1 0 sin(theta) 0 cos(theta)];

If somebody encounters similar problems, one should always observe each rotation separately and study the separate rotation matrices, for euler and axis angles, because not every function uses the same x,y,z axis or does the rotations in the standard explained way.

Anyway thanks for viewing