Problem getting slices of a volume data set along an arbitrary orientation

I have a 3D array of uint8 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 axes 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 what is 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 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

[Sean de Wolski]

You should be able to do this by applying the operations to a meshgrid of values. Generate a grid for your coordinates and then generate the z grid as a function of the x/y grid. How you get z from x/y will be based on your orientation+offset.

For example:

% Data
[x, y, z] = meshgrid(-3:0.25:3);
v = x.*exp(-x.^2 - y.^2 - z.^2);

Slice along Z = -X + Y

% Data defining a surface
[xs, ys] = meshgrid(-3:0.25:3);
zs = -xs + ys;

% Slice along it
slice(x, y, z, v, xs, ys, zs)
colorbar

[Nick]

Thank you for your answer!

The documentation is confusing for slice, it appears to only display a slice on screen, as opposed to storing it as a 2D array for me to use later. How would I go about getting this data from the slice function?

Furthermore, I am loathe to use `slice` because xs, ys and zs have to each be the size of my dataset, which is quite large when the data size is expected to be 2000x2000x2000 points.

Additionally, how should I work out xs, ys and zs, using my direction vector? I believe my current direction rotation code (and therefore the generation of my orientations and offsets) to be faulty.

[Sean de Wolski]

You're welcome, and welcome to MATLAB Answers!

First, xs/ys/zs would only have to be 2000x2000.

Second, you are correct that they just show the slice on the screen. I wasn't sure if this is what you wanted or not. However, there are two things you could do here to get the data as a matrix!

The easiest way would be to grab the CData from the slice since this is the data that is being plotted with color (i.e. the values)

s = slice(x, y, z, v, xs, ys, zs)
im = s.CData;

Alternatively, you could skip the slice command and just call griddedInterpolant to interpolate the volume for your plane. This is what slice is doing under the hood to generate the CData. http://www.mathworks.com/help/releases/R2014b/matlab/ref/griddedinterpolant-class.html

As for getting the plane from the normal, I'm going to have to think about that for a few but it shouldn't be too difficult.

[Sean de Wolski]

Actually, this looks like it does roughly what you want:

http://www.mathworks.com/matlabcentral/fileexchange/32032-extract-slice-from-volume

[Nick]

Thanks, that does indeed appear to work! I'll leave this question open for a bit just in case somebody can find where I went wrong, because I'd like to be able to use my code elsewhere, but this will work for now!