Find Pixel Coordinates of 3D Volume (DICOM) Before and After Multi-planar Reformatting (MPR) with vtkImageReslice

Support
1

Hello,

I would like to compare the value of a voxel in a DICOM volume before and after reformatting (reslicing) the volume.
I am stuck on this and hope someone may know how to do this.

Here is what I do:

  1. Read a DICOM image series using vtkDICOMImageReader (I understand the quirks of this reader, and I think I handle them appropriately thanks to the advice of @dgobbi)
  2. Perform some arithmetic on the volume.
  3. Pipe the volume through vtkImageReslice, but set the reslice axes to be the identity matrix.
    …
    slicer = vtk.vtkImageReslice()
    [ … do some flipping and casting…]
    slicer.SetInputConnection(caster.GetOutputPort())
    ident_mat, inverse_ident = vtk.vtkMatrix4x4(), vtk.vtkMatrix4x4()
    ident_mat.Identity()
    vtk.vtkMatrix4x4.Invert(ident_matrix, inverse_ident)
    slicer.SetResliceAxes(inverse_ident)
    …
  4. Display the volume, slice by slice, “in the volume plane”. That is, without reformatting. That is, without reslicing axially.

At this point, I can examine the values of the displayed pixels. Let’s say I select a pixel to examine: [x=25, y=25, z=0]. I note this pixel’s corresponding scalar value.

Next, I will apply the DICOM’s affine transformation matrix to this volume and display the reformatted (resliced) volume. (Note that the DICOM was acquired in a double-oblique orientation. So reformatting/reslicing will result in a very different-looking volume.) I want to find that pixel that I examined before.

How can I find it? Can I somehow use the DICOM’s affine transformation matrix?

Thank you in advance,
Michelle

2

What works for us every time for solving such problems, without any trial and error, is the following process:

  • Define all coordinate systems in your application: a short unique name, origin, axes unit, and axes directions. You can describe these in a table and it usually helps if you draw a sketch. See an example here. Your coordinate systems will be something like: InputVolumeIJK, InputVolumeLPS, AlignedVolumeLPS, SliceView.
  • Name all your transforms consistently, as transforms between two coordinate systems (e.g., AtoBTransform). Just by looking at your coordinate system definitions, it should be always simple to determine transformation matrix between two coordinates systems. These pair of transforms will connect all your transforms in a graph.
  • You can get transforms between any two coordinate systems by multiplying the transformation matrices that get you from one coordinate system to the other (this is where the sketch that you draw in the first step is particularly useful).