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.
The only way I found to achieve this was the Grayscale model maker module (with threshold set to 10 for the above example)
This returns a nice sphere surface, as expected, let’s call it sphere_model
Use the resulting sphere_model model and the image volume in the Probe volume with model module
this returns a model, let’s call it probed_model
The probed_model appears as a green sphere. All the surface is uniformly covered with green color.
Work the results with the models module
I tried to adjust things in the Scalars section of the module, but it has no effect on the depicted model.
What I am expecting to see is a probed_model that is colored `half-half’ with two different colors. I am new to slicer and I suppose I went wrong somewhere (and or did not yet find the right module).
I am running slicer on a MacBook (late 2013) with macOS Catalania 10.15.4 and I am using slicer version 4.10.2
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.
Here is a link with the dataset that results on my machine, following the steps I described above: probed_model.vtk - Google Drive
All the values in the scalar are set to the value of 10, which means that probably the generated volume was not right. Could you also upload the volume that you created?
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).
What you uploaded was the volume rendering property node. If you are not sure which file is which then you can click on the package icon and save the entire scene as a single .mrb file.
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.
Ok, great. I am not familiar with vtk, but it seems there are some good resources. I’ll try to learn about it and give an update once I have made some progress.
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)
def get_colored_model(index=0):
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()
contour.SetInputData(labelmapVolumeNode.GetImageData())
contour.ComputeScalarsOn()
contour.SetValue(0, index)
contour.Update()
# smoothen mesh
smoother = vtk.vtkWindowedSincPolyDataFilter()
smoother.SetInputConnection(contour.GetOutputPort())
smoother.SetNumberOfIterations(30)
smoother.NonManifoldSmoothingOn()
smoother.NormalizeCoordinatesOn()
smoother.Update()
smoothedMesh = smoother.GetOutput()
# cell normals
triangleCellNormals = vtk.vtkPolyDataNormals()
triangleCellNormals.SetInputData(smoothedMesh)
triangleCellNormals.ComputeCellNormalsOn()
triangleCellNormals.ComputePointNormalsOff()
triangleCellNormals.ConsistencyOn()
triangleCellNormals.AutoOrientNormalsOn()
triangleCellNormals.Update()
# 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)
triangleCellNormals.GetOutput().GetPointData().SetScalars(vtk_data_array)
# for mesh model
myNode = slicer.modules.models.logic().AddModel(triangleCellNormals.GetOutputPort())
# Show scalar values at surface of a model
myNode.GetDisplayNode().ScalarVisibilityOn()
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.
If you enable “Show 3D” option in the Segment Editor then you can see a 3D representation synchronized with the labelmap representation (internally quite similarly as you did it).
If you need the result as models (surface meshes) then Segment editor can create them you, just right-click on the segmentation node in Data module and choose export to model node.
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 :>