Align marching cube output to volume image using image header information

Support
1

Dear all,
I am stuck with an issue that’s probably related to the combination of VTK, SimpleITK, and 3DSlicer, any hints highly appreciated.
I am using Python.

What I am trying to do:

  1. Load a segmentation from a Nifti-File, create a mesh with vtk.vtkMarchingCubes()
  2. Do measurements on this mesh
  3. Save mesh and measurements for processing and visualization

This works fine, but when I open the vtp-files together with the original label map in 3D Slicer, they are not aligned. I am aware that Slicer uses RAS orientation, and I am most likely working with LPS. And I think, that I need to apply the affine transform from the Nifti-header to the vtkPolyData to get the correct orientation and position.
I tried the following:

  • Load the segmentation additionally with SimpleITK, get the transform matrix, apply this transform to the mesh.
  • The same as above, with inverting the first two diagonal elements of the transform matrix, to account for the LPS to RAS conversion
  • Same as above, with centering the image volume first

Any hints highly appreciated!

Best regards and thank you,
Yannick

2

If you use a very recent Slicer Preview Release then you can load a segmentation Nifti file directly as a segmentation node (choose Segmentation in the Description column in “Add data” window). Segmentation can be displayed in slice and 3D views - marching cubes (more precisely flying edges) is performed automatically and all necessary coordinate transforms are taken care of.

3

If you want, you can apply the NIFTI header transform to the mesh before you save it.

reader = vtk.vtkNIFTIImageReader()
# [set filename]

transform = vtk.vtkTransform()
if reader.GetQFormMatrix():
    transform.SetMatrix(reader.GetQFormMatrix())
elif reader.GetSFormMatrix():
    transform.SetMatrix(reader.GetSFormMatrix())

transformPoly = vtk.vtkTransformPolyData()
transformPoly.SetInputConnection(marchingcubes.GetOutputPort())
transformPoly.SetTransform(transform)
transformPoly.Update()

Note that vtkNIFTIImageReader uses RAS, so the transform described above will put the mesh into RAS.

4

Thank you for the suggestion. Unfortunately, that is exactly what I already tried. I can make it work with some examples if I flip all axes and center the volume before transforming, but I am afraid that this solution is a bit too specific and “hacky” to be widely applicable. I now reorient everything with FSL before processing and center the volumes to be on the safe side…

5

Hi Andras,
Thank you for the solution with Slicer, I am currently testing my vtk/itk approach by comparing to the output I get from Slicer. I want to process a large amount of segmentations. The calculations are quite computationally expensive, so I’m doing that on a HPC. I had a look at the code from these modules. Hope I can get to a solution that at least works for my given data.

6

If you carefully document what coordinate system is used in each data file on disk and in each data object in your code then implementing the right transformation should be straightforward.

Volume “centering” is a destructive operation (the physical position of the volume is lost), so it should be avoided, especially in the middle of a processing pipeline.

There are some ambiguities in how image orientation is defined in Nifti/Analyze formats. If you don’t specifically need to use Nifti then it might be simpler to use NRRD format, which does not have such problem.

7

Thanks for the hint with the Nifti ambiguities! I’ll look into that. With centering I did not mean the destructive Slicer function, but moving the image volume center to RAS zero and adapting the orientation matrix accordingly.
Thanks!

8

That is exactly the lossy operation I talk about. Depending on your workflow, it may not cause problems (especially if you do it before doing any other operation), but it should not be necessary. Standardization of coordinate system between various scans is an important operation but simply setting image origin to be the geometrical center of the volume extent would be a very poor solution.