Question on STL transformation in slicer's world coordinate

Support
1

Hello,

First, I created a segment and converted it into a model in slicer.

Second, I have written an external script purely based on VTK, which does not rely on Slicer classes. The script reads DICOM files using vtkDICOMImageReader, applies a threshold using vtkImageThreshold, converts it to a model using VTK’s marching cubes, and finally exports it using vtkSTLWriter.

Now, I imported the STL generated by the script, into Slicer to compare it against the model generated earlier from the segment. However, I am unable to compare their contours in 2D because my external STL is not aligned with the model generated by Slicer. In the 3D scene, both models are far apart and do not lie on the same origin. I don’t want to use the Transform module because of its manual nature. I want to simulate, how slicer internally does transformations, within my script so that when I export any STL from my script and import it into Slicer, it has the correct transformations and aligns with the 3D Slicer-generated model.

2

https://slicer.readthedocs.io/en/latest/user_guide/coordinate_systems.html

3

Thank you for providing the reference.

I have gone through it, and I was wondering how I can utilize the exact transformation matrix of Slicer. I am setting the labelmap’s origin and spacing to be the same as the one given in the CreateClosedSurface function of the vtkBinaryLabelmapToClosedSurfaceConversionRule class.

I’ve been stuck on this problem for many days. I would appreciate further help.

Here’s a picture of the models for reference (the green model is from Slicer, and the yellow one is generated from the script):

image
image326×246 8.46 KB

Thank you.

4

My guess is that you just need to pick LPS on loading like this:

image
image930×324 30.2 KB

If that’s it and you want to make this automatic in the future you can add a coordinate system tag in the stl header something like this:

https://slicer.readthedocs.io/en/v4.11/user_guide/data_loading_and_saving.html#models

5

Hello @pieper ,

I tried loading the STL by selecting both coordinate systems, but it didn’t work, as the result is still the same. I also modified my script to include the coordinate tag, and when I loaded the STL again, there was no luck.

Here’s my code for saving a mesh:

def SaveMeshToSTL(polydata, filename):
    stl_writer = vtk.vtkSTLWriter()
    stl_writer.SetFileName(filename)

    stl_writer.SetInputData(polydata)

    coordinateSystemTag = "SPACE"
    CoordinateSystem = slicer.vtkMRMLModelStorageNode().CoordinateSystemLPS
    coordinateSystemStr = slicer.vtkMRMLModelStorageNode().GetCoordinateSystemTypeAsString(CoordinateSystem)
    coordinateSytemSpecification = coordinateSystemTag + "=" + coordinateSystemStr
    header = "3D Slicer output. " + coordinateSytemSpecification
    stl_writer.SetHeader(header)

    stl_writer.SetFileType(vtk.VTK_BINARY)

    stl_writer.Write()
    print(f"Model successfully saved to {filename}")

The output of the function TransformIJKToWorld is then passed to the mesh-saving function as polydata:

def TransformIJKToWorld(polydata, oriented_data):
    labelmapGeometryTransform = vtk.vtkTransform()
    labelmapImageToWorldMatrix = vtk.vtkMatrix4x4()

    oriented_data.GetImageToWorldMatrix(labelmapImageToWorldMatrix)
    labelmapGeometryTransform.SetMatrix(labelmapImageToWorldMatrix)

    transformPolyDataFilter = vtk.vtkTransformPolyDataFilter()
    transformPolyDataFilter.SetInputData(polydata)
    transformPolyDataFilter.SetTransform(labelmapGeometryTransform)
    transformPolyDataFilter.Update()

    return transformPolyDataFilter.GetOutput()

@lassoan @cpinter @muratmaga, Any idea what might be the issue?

6

Well, it must be something else then. You just need to trace through all the matrices and coordinates to validate your assumptions.

1 Like
7

vtkDICOMImageReader only works for special cases. I would not recommend spending any time with trying to debug why the image is not loaded correctly with vtkDICOMImageReader. If you are looking for a pure VTK-based solution then your best option is GitHub - dgobbi/vtk-dicom: A set of classes for using DICOM in VTK.. It is still very limited but works for majority of basic image acquisitions. Make sure you take into account the PatientMatrix provided by the reader.

If you want to reconstruct volumes from more complex acquisitions then you can use dcm2niix or Slicer. If you want an open-source solution to read not just images but other important DICOM objects then probably your only choice is Slicer (or implement yourself from pydicom level).

1 Like
8

@lassoan, your suggested solution really worked! I used the PatientMatrix from vtkDICOMReader, and it’s giving me correct results now.

Thank you!