Calculation of 4x4 matrix of planes

Hi,
I am using the example of ImplicitPlaneWidget2.
I want to be able to manipulate this by inputting a 3D image (.raw).

I studied about 4x4 matrix.
I would like to write a program to calculate 4x4 matrix in plane.

The code in vtkImageResliceMapper can not understand the following “axis, saxis, taxis”.
Can you explain it?

// Create two orthogonal axes
double saxis[3], taxis[3];
taxis[0] = 0.0;
taxis[1] = 1.0;
taxis[2] = 0.0;
if (maxi == 1)
{
  taxis[1] = 0.0;
  taxis[2] = 1.0;
}
vtkMath::Cross(taxis, axis, saxis);

https://gitlab.kitware.com/vtk/vtk/blob/master/Rendering/Image/vtkImageResliceMapper.cxx#L600

@dgobbi

It is hard to explain without pictures, but I will try.

The vectors saxis[3], taxis[3], and axis[3] are three orthogonal axes where ‘axis’ is roughly aligned with the slice normal. The only important thing about ‘saxis’ and ‘taxis’ is that they are orthogonal to ‘axis’.

The code below computes the rotation between ‘axis’ and ‘normal’:

// Compute the rotation angle between the axis and the normal
double vec[3];  // vec[3] is the axis of rotation
vtkMath::Cross(axis, normal, vec);
double costheta = vtkMath::Dot(axis, normal);
double sintheta = vtkMath::Norm(vec);
double theta = atan2(sintheta, costheta);

Then ‘vec[3]’ and ‘angle’ are used to compute a 3x3 rotation matrix. We use this matrix to rotate ‘saxis[3]’ and ‘taxis[3]’:

// Create a slice-to-data transform matrix
// The columns are v1, v2, normal
double v1[3], v2[3];
vtkMath::Multiply3x3(mat, saxis, v1);
vtkMath::Multiply3x3(mat, taxis, v2);

resliceMatrix->Element[0][0] = v1[0];
resliceMatrix->Element[1][0] = v1[1];
resliceMatrix->Element[2][0] = v1[2];
resliceMatrix->Element[3][0] = 0.0;

resliceMatrix->Element[0][1] = v2[0];
resliceMatrix->Element[1][1] = v2[1];
resliceMatrix->Element[2][1] = v2[2];
resliceMatrix->Element[3][1] = 0.0;

resliceMatrix->Element[0][2] = normal[0];
resliceMatrix->Element[1][2] = normal[1];
resliceMatrix->Element[2][2] = normal[2];
resliceMatrix->Element[3][2] = 0.0;

The whole purpose of this is to compute vectors ‘v1[3]’ and ‘v2[3]’ that are perpendicular to ‘normal[3]’. All of this fancy math is used to ensure that ‘v1[3]’ and ‘v2[3]’ are as close as possible to the original volume axes.

Thank you for your reply.
Sorry for my late reply.

I could understand about 4x4 matrix by your detailed explanation.

May I ask another question?
I’m having a hard time understanding the linked code from line 600 to line 642.
Could you explain briefly?

@dgobbi

The ‘plane[4]’ is the coefficients a,b,c,d for the
equation of a plane. It is important to understand that the first three components a,b,c of plane[4] are the normal, so on line 664 it is possible to say normal = plane.

LInes 600 to 630 get the plane (and wplane, which is the plane after the actor transformation from ‘data’ coords to ‘world’ coords). The plane is oriented so that the normal points towards the ‘camera’ (this means it points out of the computer screen, not into the computer screen).

Lines 631 to 642 do exactly what the comment says: find the largest component of the normal. This checks whether the normal is closest to the x, the y, or the z axis.

For your original question, I think it is most important to understand how ‘v1’ and ‘v2’ are computed. Compare this to the simpler (much less general) answer I gave to Post 1391, which computes the matrix only for a rotation around the Z axis:

Thank you for your detailed explanation.
It helped me to understand.

Thanks to your help, I have a rough understanding of the cut plane matrix and compute method.
Also, I could understand how to compute ‘v1’ and ‘v2’.

What should I do as a next step in creating a program?
Is it difficult to incorporate it into the program of ImplicitPlaneWidget2 from the beginning?

I’m having trouble getting started with program improvement.

@dgobbi

You’re asking a little too much. Just start anywhere, and when you start things will become clearer. When you have some code written, I’d be glad to check it for you.

Thank you for your reply.

I will try to make a program by myself.
At that time, I would appreciate your advice.

Based on the code of vtkImageResliceMapper.cxx, I created a program to calculate 4 x 4 matrix on a plane.

The code created is shown below.

  // Setup a visualization pipeline
  vtkSmartPointer<vtkPlane> plane =vtkSmartPointer<vtkPlane>::New();
  double normal[3];
  plane->GetNormal(normal);
  double Origin[3];
  plane->GetOrigin(Origin);
  // Set the slice orientation
 vtkSmartPointer<vtkMatrix4x4> resliceMatrix =
vtkSmartPointer<vtkMatrix4x4>::New();

// Create the corresponding axis
double axis[3];
axis[0] = 0.0;
axis[1] = 0.0;
axis[2] = 0.0;

double saxis[3], taxis[3];
taxis[0] = 0.0;
taxis[1] = 1.0;
taxis[2] = 0.0;
vtkMath::Cross(taxis, axis, saxis);

// Compute the rotation angle between the axis and the normal
double vec[3];
vtkMath::Cross(axis, normal, vec);
double costheta = vtkMath::Dot(axis, normal);
double sintheta = vtkMath::Norm(vec);
double theta = atan2(sintheta, costheta);
if (sintheta != 0)
{
  vec[0] /= sintheta;
  vec[1] /= sintheta;
  vec[2] /= sintheta;
}

// convert to quaternion
costheta = cos(0.5*theta);
sintheta = sin(0.5*theta);
double quat[4];
quat[0] = costheta;
quat[1] = vec[0]*sintheta;
quat[2] = vec[1]*sintheta;
quat[3] = vec[2]*sintheta;

// convert to matrix
double mat[3][3];
vtkMath::QuaternionToMatrix3x3(quat, mat);

// Create a slice-to-data transform matrix
// The columns are v1, v2, normal
double v1[3], v2[3];
vtkMath::Multiply3x3(mat, saxis, v1);
vtkMath::Multiply3x3(mat, taxis, v2);

resliceMatrix->Element[0][0] = v1[0];
resliceMatrix->Element[1][0] = v1[1];
resliceMatrix->Element[2][0] = v1[2];
resliceMatrix->Element[3][0] = 0.0;

resliceMatrix->Element[0][1] = v2[0];
resliceMatrix->Element[1][1] = v2[1];
resliceMatrix->Element[2][1] = v2[2];
resliceMatrix->Element[3][1] = 0.0;

resliceMatrix->Element[0][2] = normal[0];
resliceMatrix->Element[1][2] = normal[1];
resliceMatrix->Element[2][2] = normal[2];
resliceMatrix->Element[3][2] = 0.0;

resliceMatrix->Element[0][3] = Origin[0];
resliceMatrix->Element[1][3] = Origin[1];
resliceMatrix->Element[2][3] = Origin[2];
resliceMatrix->Element[3][3] = 1.0;

Can you tell me if there is something to be improved?

@dgobbi

You misunderstood the purpose of ‘axis’, ‘taxis’ and ‘saxis’. Here are two concrete examples that might help. The goals are:

  • ‘axis’ must be approximately the same direction as ‘normal’
  • ‘saxis’, ‘taxis’, and ‘axis’ must be orthonormal (and must obey right-hand-rule)

Example 1) if ‘normal’ is [0.48716132, 0.84723708, -0.21180927] then:

  • ‘axis’ should be [0, 1, 0]
  • ‘taxis’ should be [0, 0, 1]
  • ‘saxis’ should be [1, 0, 0]

Example 2) if ‘normal’ is [0.2470483, 0.19003716, -0.95018578] then:

  • ‘axis’ should be [0, 0,-1]
  • ‘taxis’ should be [0, 1, 0]
  • ‘saxis’ should be [1, 0, 0]

The purpose of the rest of the math is to constrain the rotation of the 4x4 matrix such that the axis-of-rotation is parallel to the plane (so there is no in-plane rotation).

Thank you for your reply.

I could understand ‘axis’, ‘taxis’ and ‘saxis’ by example.
It has been improved as follows.

int maxi = 0;
double maxv = 0.0;
for (int i = 0; i < 3; i++)
{
  double tmp = normal[i]*normal[i];
  if (tmp > maxv)
  {
    maxi = i;
    maxv = tmp;
  }
}
// Create the corresponding axis
double axis[3];
axis[0] = 0.0;
axis[1] = 0.0;
axis[2] = 0.0;
axis[maxi] = ((normal[maxi] < 0.0) ? -1.0 : 1.0);
double saxis[3], taxis[3];
taxis[0] = 0.0;
taxis[1] = 1.0;
taxis[2] = 0.0;
if (maxi == 1)
{
  taxis[1] = 0.0;
  taxis[2] = 1.0;
}
vtkMath::Cross(taxis, axis, saxis);

Is the fourth row of the 4x4 matrix correct?

resliceMatrix->Element[0][3] = Origin[0];
resliceMatrix->Element[1][3] = Origin[1];
resliceMatrix->Element[2][3] = Origin[2];
resliceMatrix->Element[3][3] = 1.0;

Before I answer that question, I will give a brief summary. So far, we have computed a rotation that takes the 2D plane (the XY plane at Z=0, which will be the output of vtkImageReslice) and rotates it so that it becomes parallel with the 3D slice plane. The next thing we have to do is translate the rotated 2D plane so that is in the same position as the 3D slice plane.

The fourth column of the 4x4 matrix gives the translation.

So the quick answer to your question is: yes, setting the translation to the plane “Origin” is a solution. This works because the point (0,0,0) is on the 2D plane (obviously!), and if we rotate that point it is still (0,0,0), and if we then translate that point by (Origin[0], Origin[1], Origin[2]), then it moves to the vtkPlane’s Origin.

A longer answer is that setting the translation to “Origin” is a bad solution. It is bad because (0,0,0) is usually the corner of the 2D plane (the vtkImageReslice output), and “Origin” is usually the center of the slice plane (controlled by vtkImplicitPlaneWidget2).

It is better to translate the center of the 2D plane to the center of the slice plane.

Let’s define the following variables:

  • q = center of 2D plane (center computed from the dimensions of the slice)
  • p = center of 3D slice plane (Origin of the vtkPlane)
  • R = 3x3 rotation matrix that we computed previously
  • t = translation

And write this transformation equation:

  • p = Rq + t

Solve for t:

  • t = p - Rq

This gives the general solution, if we expand the matrix multiplication:

resliceMatrix->Element[0][3] = origin[0] - center[0]*v1[0] - center[1]*v2[0] - center[2]*normal[0];
resliceMatrix->Element[1][3] = origin[1] - center[0]*v1[1] - center[1]*v2[1] - center[2]*normal[1];
resliceMatrix->Element[2][3] = origin[2] - center[0]*v1[2] - center[1]*v2[2] - center[2]*normal[2];
resliceMatrix->Element[3][3] = 1.0;

In case you are wondering how to compute ‘center’, the answer is that first you must decide on the dimensions of your output slice. For example, if your input volume is 512x512 and voxel spacing is 1 millimeter, then you will probably want to use something like this:

int n = 512; // size in pixels
double s = 1.0; // size of each pixel
reslice->SetOutputOrigin(0.0, 0.0, 0.0);
reslice->SetOutputExtent(0, n-1, 0, n-1, 0, 0);
reslice->SetOutputSpacing(s, s, s);

double center[3] = { 0.5*(n-1)*s, 0.5*(n-1)*s, 0.0 };

You can choose ‘n’ to be whatever you want. If you want to avoid cropping any of the volume when you slice it, then make ‘n’ to be 1.5 times the largest volume dimension.

Hi David,
Thank you for your advice.

I was able to understand the translation of the plane.

Should I apply a matrix to the following vtkImageReslice after this?

vtkSmartPointer<vtkImageReslice> reslice =
vtkSmartPointer<vtkImageReslice>::New();
reslice->SetInputConnection(reader->GetOutputPort());
reslice->SetOutputDimensionality(2);
reslice->SetResliceAxes(resliceMatrix);
reslice->SetInterpolationModeToLinear();
reslice->SetOutputOrigin(0.0, 0.0, 0.0);
reslice->SetOutputExtent(0, n-1, 0, n-1, 0, 0);
reslice->SetOutputSpacing(s, s, s);

Also in the last comment, could you explain a little more about about 1.5 times the maximum volume?

@dgobbi

Maximum volume dimension. So for a 200x500x300 volume, the maximum dimension is 500 (assume isotropic spacing).

The question is this: what are the maximum possible dimensions for an oblique slice, if the oblique rotation is unknown? Obviously larger than 200x500 or 500x300. Even 500x500 might not be large enough to contain the slice. But 750x750 will almost definitely be large enough.

The answer is that the dimensions of an oblique slice can be (roughly) up to 1.5 times the dimensions of the original slices. Usually this is an overestimate.

Thank you for your explane.
Sorry for my late reply.

I was able to understand my question with the detailed explanation.