DICOM MRI Multiplanar Reformation (MPR) using vtkImageReslice

Hello,

I am new to VTK. I’m stuck on a problem and could really use some advice. I am writing a small program to perform multiplanar reformation (MPR) on a DICOM. I’ve read many forums, Q&As (@dgobbi) and examples. And my program is almost there, but not quite.

Here is what I am trying to do:
Given two DICOM-formatted MRIs - one which was acquired orthogonally and the other acquired in any orientation (usually doubly oblique), reformat the oblique MRI to display overlain on the orthogonal MRI. Note that the two volumes will likely have different resolutions (different voxel sizes). Also, the two volumes will always overlap. That is, in patient space, the oblique MRI lies within the volume of the orthogonal MRI.

Here is how I am trying to do it:

  1. Use vtkDICOMImageReader to load two DICOM MRI series. One MRI acquired orthogonally, and the other acquired in a doubly oblique orientation.
  2. Use ImageOrientationPatient and ImagePositionPatient from each DICOM to create an affine transformation matrix that describes the orientation of each MRI in patient space.
  3. Create a vtkImageReslice object for each MRI, and use the affine transformation matrices to set the reslice axes.
  4. Display both reformatted MRI’s, presumably now overlapping, using vtkRenderer. I display the volumes one slice at a time, with a mouse-move callback to allow scrolling through slices.

The problem:
When the resulting two volumes display in the render window, the reformatted oblique volume is in the wrong orientation, and displays far outside the orthogonal volume.

What am I doing wrong? Do I need to be somehow using SetOutputOrigin() with the reslice objects?
Thank you so much to anyone who reads all the way through this post!

  • Michelle K.

My Code:

Summary

import vtk

===== IMAGE 1 (orthogonal) =====

image1_reader = vtk.vtkDICOMImageReader()
[…set directory name, update reader…]

vtkDICOMImageReader sorts images in reverse order. Flip front-to-back.

image1_flip = vtk.vtkImageReslice()
image1_flip.SetInputConnection(image1_reader.GetOutputPort())
image1_flip.SetResliceAxesDirectionCosines(1, 0, 0, 0, 1, 0, 0, 0, -1)
image1_flip.Update()

create affine matrix

image1_position = image1_reader.GetImagePositionPatient()
image1_orient = image1_reader.GetImageOrientationPatient()
image1_orient_x = image1_orient[0:3]
image1_orient_y = image1_orient[3:6]
math = vtk.vtkMath()
image1_orient_z = [0 for i in range(3)]
math.Cross(image1_orient_x, image1_orient_y, image1_orient_z)
image1_affine = vtk.vtkMatrix4x4()
for i in range(0, 3):
image1_affine.SetElement(i, 0, image1_orient_x[i])
image1_affine.SetElement(i, 1, image1_orient_y[i])
image1_affine.SetElement(i, 2, image1_orient_z[i])
image1_affine.SetElement(i, 3, position[i])

create a reslice object for the image 1

image1_reslice = vtk.vtkImageReslice()
image1_reslice.SetInputConnection(image1_flip.GetOutputPort())

Set the slicing axis for image 1

image1_reslice.SetOutputDimensionality(2)
image1_reslice.SetResliceAxes(image1_affine)
image1_reslice.SetInterpolationModeToLinear()

===== IMAGE 2 =====

[… repeat for image 2 …]

===== Map and Render Both Images =====

renderer = vtk.vtkRenderer()
image1_mapper = vtk.vtkImageResliceMapper()
image1_mapper.SetInputConnection(image1_reslice.GetOutputPort())
image1 = vtk.vtkImageSlice()
image1.SetMapper(image1_mapper)
renderer.AddViewProp(image1)

image2_mapper = vtk.vtkImageResliceMapper()
image2_mapper.SetInputConnection(image2_reslice.GetOutputPort())
image2 = vtk.vtkImageSlice()
image2.SetMapper(image2_mapper)
image2.GetProperty().SetOpacity(0.7)
renderer.AddViewProp(image2)

renWin = vtk.vtkRenderWindow()
renWin.SetSize(750, 750)
renWin.AddRenderer(renderer)
style = vtk.vtkInteractorStyle()
iren = vtk.vtkRenderWindowInteractor()
iren.SetInteractorStyle(style)
renWin.SetInteractor(iren)
renWin.Render()

Create callbacks for slicing the image

actions = {}
actions[“Slicing”] = 0

def button_callback(obj, event):

def mouse_move_callback(obj, event):

style.AddObserver(“MouseMoveEvent”, mouse_move_callback)
style.AddObserver(“LeftButtonPressEvent”, button_callback)
style.AddObserver(“LeftButtonReleaseEvent”, button_callback)

Start interaction

iren.Start()

Hi Michelle, welcome to the board.

Here is some old code that I wrote for reading an oriented image with vtkDICOMImageReader (the whole program is here). This function expects the user to create an empty image and matrix (with the usual ::New()` method), and will fill them with the DICOM slices and orientation matrix.

In particular, note that the code sets the Origin to (0,0,0) and stores the position in the last column of the matrix.

void ReadDICOMImage(
  vtkImageData *data, vtkMatrix4x4 *matrix, const char *directoryName)
{
  // read the image
  vtkSmartPointer<vtkDICOMImageReader> reader =
    vtkSmartPointer<vtkDICOMImageReader>::New();

  reader->SetDirectoryName(directoryName);
  reader->Update();

  // the reader flips the image and reverses the ordering, so undo these
  vtkSmartPointer<vtkImageReslice> flip =
    vtkSmartPointer<vtkImageReslice>::New();

  flip->SetInputConnection(reader->GetOutputPort());
  flip->SetResliceAxesDirectionCosines(
    1,0,0, 0,-1,0, 0,0,-1);
  flip->Update();

  vtkImageData *image = flip->GetOutput();

  // get the data
  data->CopyStructure(image);
  data->GetPointData()->PassData(image->GetPointData());
  data->SetOrigin(0,0,0);

  // generate the matrix
  float *position = reader->GetImagePositionPatient();
  float *orientation = reader->GetImageOrientationPatient();
  float *xdir = &orientation[0];
  float *ydir = &orientation[3];
  float zdir[3];
  vtkMath::Cross(xdir, ydir, zdir);

  for (int i = 0; i < 3; i++)
    {
    matrix->Element[i][0] = xdir[i];
    matrix->Element[i][1] = ydir[i];
    matrix->Element[i][2] = zdir[i];
    matrix->Element[i][3] = position[i];
    }
  matrix->Element[3][0] = 0;
  matrix->Element[3][1] = 0;
  matrix->Element[3][2] = 0;
  matrix->Element[3][3] = 1;
  matrix->Modified();
}

I definitely consider the above to be “old” code. Currently I use vtkDICOMReader from vtk-dicom instead, since it provides a much more complete DICOM implementation and is has a method SetMemoryRowOrderToFileNative() that ensures proper orientation in LPS coords.

When displaying the images, it’s best to use vtkImageStack to ensure that VTK does the image compositing properly. The vtkImageSlice objects go into the stack, and then the stack goes into the renderer. Here is an example: https://lorensen.github.io/VTKExamples/site/Cxx/Images/ImageStack/

Usually I do not use vtkImageReslice to handle the orientation, since the orientation can be handled automatically by vtkImageResliceMapper if pass the matrix to the vtkImageSlice, e.g. image1->SetUserMatrix(image1_affine). You can take a look at the program that I referred to above for an example of how to do this.

If you do decide apply the affine transformation with vtkImageReslice instead of with the mapper/actor via SetUserMatrix(), then please note that the ResliceAxes matrix specifies a transformation from the reslice output coords to reslice input coords. So in your sample code, you would have to invert the matrix. This is another reason why SetUserMatrix() can be an easier way to supply the affine matrix, since it takes the forward transformation rather than the reverse transformation.

In situation where I do want to use vtkImageReslice to handle the orientation for an overlay, then what I usually do is reslice my ‘overlay’ image to match my ‘base’ image, rather than reslicing both of the images.

1 Like

Hi David,

Thank you for your thoughtful answer. I truly appreciate all of your responses to vtk users’ questions.

Since I am a new vtk user, it took me a little bit to thoroughly understand your suggestions. My script works now. The key solution for me was realizing that the y- and z- axis flipping that the vtkDICOMImageReader quietly does is significant. If I do the following:

  1. read the DICOM
  2. flip image in Y
  3. flip the image in Z
  4. construct the affine
  5. inverse the affine
  6. pipe the (y- and z- flipped) image into a vtkImageReslice object
  7. apply the affine to the vtkImageReslice object using SetResliceAxes
    … then the resulting image is properly aligned (orthogonal / axial).

I don’t have access to the vtkDICOMReader package, because I am running the vtk libraries from a downloaded executable. I did not build vtk from source code.

I am going to try using vtkImageResliceMapper and also the vtkImageStack, as you recommend.

Michelle