Get slices of a volume data set along an arbitrary orientation

Question

I have a 3D array of integer values, and a direction vector. I want to be able to pick out X slices along this vector, with the orientation being the direction the slicer is looking (that is to say, the X and Y coordinates of the slice are two normals to the direction vector that are perpendicular to each other).

I have the following solution, but the results show at a 90-degree offset to expected along the Y axis local to the direction vector, and I cannot work out why.

%Get "up" and "right" direction vectors
up = get_rotated_direction(direction, 90,0,0);
right = get_rotated_direction(direction, 0,90,0);

z = 0;

for i=1:zStep:zCount
    z = z +1;
    %Slicer follows a path denoted by "points"
    point = points(i,:);

    %topLeft holds the coordinates in the data volume of the point at [X=0,Y=0] in the slice.
    topLeft = point - right * width / 2 - up * height / 2;

    for x=1:width
        for y=1:height
            offsetPoint = round(topLeft + up*y+ right*x);
            result(x,y,z) = data(offsetPoint(1),offsetPoint(2),offsetPoint(3));
        end
    end
end

The result is a valid image, but at the wrong orientation, so I believe the issue lies in get_rotated_direction or my calls to it, however in testing with the unit vectors [0,0,1], [0,1,0] and [1,0,0] the results appear to be correct.

get_rotated_direction looks like this:

function result = get_rotated_direction (direction, x_angle,y_angle,z_angle)
    %get the rotation from origin to direction
    rot = vrrotvec([0,0,1],direction);

    %Turn this into a quaternion
    rotquat = angle2quat(rot(1)*rot(4),rot(2)*rot(4),rot(3)*rot(4),'XYZ');

    %get the angles in radians
    x_angle = deg2rad(x_angle);
    y_angle = deg2rad(y_angle);
    z_angle = deg2rad(z_angle);

    %Create the desired rotation from the origin
    desiredRot = angle2quat(x_angle,y_angle,z_angle,'XYZ');

    %Rotate the origin to our desired rotation
    %Then rotate that by the rotation from origin to our direction
    %This should produce a direction that has been locally rotated.
    direction= quatrotate(rotquat,quatrotate(desiredRot,[0,0,1]));

    %Normalize the direction and return
    result = direction/norm(direction);
end

Answer 1

Theres a more general approach for this using transformation matrices. I had a similar problem and used this answer, which I tried to generalise to the 3D case below:

%example data
load mri; olddata = squeeze(D);
phi = 90; theta = 90; psi = 0;
r = [phi, theta, psi]; r = r*(pi/180);
dims = size(olddata);
p0 = [round(dims(1)/2);round(dims(2)/2);round(dims(3)/2)]; %image center

%I assume that spacing between points is equal in all three dimensions
Rx=[1 0 0 0;
    0 cos(r(1)) sin(r(1)) 0;
    0 -sin(r(1)) cos(r(1)) 0;
    0 0 0 1];
Ry=[cos(r(2)) 0 -sin(r(2)) 0;
    0 1 0 0;
    sin(r(2)) 0 cos(r(2)) 0;
    0 0 0 1];
Rz=[cos(r(3)) sin(r(3)) 0 0;
    -sin(r(3)) cos(r(3)) 0 0;
    0 0 1 0;
    0 0 0 1];
R = Rz*Ry*Rx;
%make affine matrix to rotate about center of image instead of origin
T = ( eye(3)-R(1:3,1:3) ) * p0(1:3);
A = R;
A(1:3,4) = T;
Rold2new = A;
Rnew2old = inv(Rold2new);

%this is the transformation
%I assume you want a volume with the same dimensions as your old volume
[xx yy zz] = meshgrid(1:dims(1),1:dims(2),1:dims(3));
coordinates_axes_new = [xx(:)';yy(:)';zz(:)'; ones(size(zz(:)))']; % the coordinates you want to find a value for
coordinates_axes_old = Rnew2old*coordinates_axes_new; % the coordinates where you can find those values
Xcoordinates = reshape(coordinates_axes_old(1,:), dims(1), dims(2), dims(3));
Ycoordinates = reshape(coordinates_axes_old(2,:), dims(1), dims(2), dims(3));
Zcoordinates = reshape(coordinates_axes_old(3,:), dims(1), dims(2), dims(3));

%this is where you get the values for your new coordinates
method = 'nearest'; %you have integers only so you probably want the nearest neighbour?
newdata = interp3(olddata, Xcoordinates, Ycoordinates, Zcoordinates, method);
figure; subplot(1,2,1); imshow(olddata(:,:,13), []); title('original'); subplot(1,2,2); imshow(newdata(:,:,13), []); title('rotated')

Which results in this rotation: