I have written a couple of python scripts to analyze CT scan images and would like to depict the results in slicer.
My script returns a 3d array that assigns a thickness to each point on a surface (let’s say the thickness is represented as an integer in [0,255], if that is important). The first thing I want to achieve is to depict the colored surface.
For the sake of simplicity let’s assume the surface is a sphere and one half of it should be colored in one color and the other in another color, according to some color-mapping. Here is an example that creates the array.
For trivial segmentation tasks (e.g., bone segmentation from CT), simple thresholding works, so “Grayscale model maker” module is usable. For most other cases, you can use Segment Editor module.
Probing with model adds a new scalar array to a model. You can use Models module’s Scalars section to make an array colorize the surface. If you have trouble setting this up then upload an example model somewhere (the result of probing) and post the link here.
Here’s an example that may help. The purpose is to show a way to use C++ and python together, but it also illustrates all the mechanism needed to manipulate model scalar arrays from numpy arrays and manage the color mapping.
Yes, indeed. The segmentation is not the problem I think. The problem is the coloring of the resulting surfaces. The sphere example I provided is just an instructive example to reproduce the issues I am facing with a simpler example.
This looks indeed helpful for what I have in mind. However, for the moment seems a bit overkill in order to understand the issue at hand for me. So I have tried to take a step back and just extract the scalar values. The result is insightful in some sense.
Not, this is all looks good. Mesh processing modules may or may not add surface normals to the output. If you need surface normals and the modules that you used did not happen to provide them then you can apply vtkPolyDataNormals to the output (or if you want to use GUI, then use Surface Toolbox module).
The scene file shows that you probed the same volume that you contoured! I have done this a few times before but of course this can never work, as per definition, the isosurface that you extract from a volume will have the same scalar value at all the surface points (the exact value that you specified for isosurface extraction). Since you have created your model by extracting the iso-value of 10, all the probed point scalars have the value of 10, which means both iso-surface extraction and probing worked perfectly.
See this example how to create a surface mesh using a VTK source. You can use a vtkSphereSource to create a sphere model.
All right, thanks! I thought that by applying the Probe volume with model module, the surface values would change according to the values of the volume. But I clearly misunderstood.
Then back to the original problem: how can I achieve what I want to do in slicer? If we stick to the example I gave. How do I extract the surface model and then color this surface model according to the entries in the array that defines the surface?
You carefuly built your mesh so that each point lies in positions that all have the exact same scalar value. So, probing will always give you that same scalar value for each point.
You did it well. Just don’t create your model by iso-surfacing the same volume that you use for coloring.
Extracting a medial surface from a volume is not a trivial operation. If during your skeletonization operation you have not built a data structure that preserved your original topology then you need to reconstruct the surface from point cloud in the end. VTK has many tools for this. For example, you can convert voxel data to points, remove zero-valued points using thresholding, then create a surface using one of the surface reconstruction filters. These all should not be more than 10-20 lines of Python code in total.
I have had my try on some vtk. I manage to do what I intended to do in vtk now, but I am failing to incorporate it properly in slicer.
Here is a minimalistic example of my code that should run from within slicer (after installing scipy, which I did as described in the forum):
import numpy as np
from scipy.ndimage.filters import maximum_filter
from vtk.util.numpy_support import vtk_to_numpy
from vtk.util.numpy_support import numpy_to_vtk
def ball(radius, dtype=np.uint8):
n = 2 * radius + 1
Z, Y, X = np.mgrid[-radius:radius:n * 1j,
-radius:radius:n * 1j,
-radius:radius:n * 1j]
s = X ** 2 + Y ** 2 + Z ** 2
return np.array(s <= radius * radius, dtype=dtype)
segmentationNode = slicer.mrmlScene.GetFirstNodeByClass("vtkMRMLSegmentationNode")
referenceVolumeNode = slicer.mrmlScene.GetFirstNodeByClass("vtkMRMLScalarVolumeNode")
labelmapVolumeNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLLabelMapVolumeNode')
# extract mask
slicer.modules.segmentations.logic().ExportVisibleSegmentsToLabelmapNode(segmentationNode, labelmapVolumeNode, referenceVolumeNode)
# apply marching cubes
contour = vtk.vtkDiscreteMarchingCubes()
# smoothen mesh
smoother = vtk.vtkWindowedSincPolyDataFilter()
smoothedMesh = smoother.GetOutput()
# cell normals
triangleCellNormals = vtk.vtkPolyDataNormals()
# get mask from segmentation
labelmapArray = slicer.util.arrayFromVolume(labelmapVolumeNode)
# extract volume
referenceVolumeArray = slicer.util.arrayFromVolume(referenceVolumeNode)
# mask original volume
# maskedVolumeArray = np.copy(labelmapArray*referenceVolumeArray)
maskedVolumeArray = labelmapArray * referenceVolumeArray
maskedVolumeArray = maskedVolumeArray.transpose(2,1,0)
# map volume values to surface
coords = vtk_to_numpy(triangleCellNormals.GetOutput().GetPoints().GetData())
voi_array = maximum_filter(maskedVolumeArray, footprint=ball(5))
mn = voi_array[voi_array > 0].min()
mx = voi_array[voi_array > 0].max()
voi_array = (voi_array-mn) / (mx-mn)
vals = voi_array[tuple(coords.astype(np.int).T)]
vtk_data_array = numpy_to_vtk(vals)
# for mesh model
myNode = slicer.modules.models.logic().AddModel(triangleCellNormals.GetOutputPort())
# Show scalar values at surface of a model
And here is the example data:
The workflow would be the following:
load example.nrrd in slicer
open the segment editor and extract the segment by applying the threshold operation.
run the get_colored_model function from the python interactor within slicer
The result is a correctly colored mesh as I intended
I am doing things without really understanding the slicer api and I am not very happy with my solution for the following reasons:
a) The model appears as an independent model, not linked to the original volume/segment (in the 3D viewer it appears as a kind of mirror-version, next to the original volume/segment).
For the moment this is fine, but it would be nicer to have the model be linked to the rest. I am planing to write a small extension to slicer that can automatically apply the algorithm.
b) I cannot interact with the colors of the resulting model in the Models module. When I try to change the color table or anything else, the depiction breaks and suddenly depicts nonsense. Ideally I would like to be able to change the color table from within slicer.
Is it possible to adapt the code to incorporate the last points?
Thanks in advance!
I resolved (b). There was a problem with the color mapping and also a bug because I got too used to Python3 (however it did not play a role for the specific example above). voi_array = (voi_array-mn) / (mx-mn) should be replaced by voi_array = (voi_array-mn) / float(mx-mn) * 255 to resolve the coloring issues.
I could not find anything in the script repository to help. I suppose I am doing sth wrong somewhere.
Edit: I have also tried to apply again the Probe Volume With Model module, but it does not resolve the issue. It gives me a proper array of coordinates, but since they are in a different coordinate system it still does not allow to apply my mapping for the Scalar Values.
Sorry, I am a bit lost in between different slicer-objects :>