# Coloring surfaces in slicer

**URL:** <https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202>\
**Category:** Support\
**Tags:** color\
**Created:** [April 20, 2020, 1:47pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202 "2020-04-20T13:47:37Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![hcmidt](https://avatars.discourse-cdn.com/v4/letter/h/59ef9b/32.png) [@hcmidt](https://discourse.slicer.org/u/hcmidt)\
**Post date:** [April 20, 2020, 1:47pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/1 "2020-04-20T13:47:37Z")

</div>

Hello together,

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.

```auto
import numpy as np
from skimage.morphology import ball

b1 = ball(100)
b2 = ball(99)
b = np.copy(b1)
b[1:-1, 1:-1, 1:-1] -= b2
b[0:100][b[0:100] != 0] = 100
b[100:][b[100:] != 0] = 250

```

I want the half that is assigned the values `100` as one color and the other half that is assigned the values `250` in another color.

For now I am just storing this data in nrrd format and load it into slicer.

```auto
import itk
image = itk.GetImageFromArray(b)
itk.imwrite(image, '~/image.nrrd')

```

I have _tried_ to use the `Probe volume with model` module (like suggested in [How to analyze the thickness of the model](https://discourse.slicer.org/t/how-to-analyze-the-thickness-of-the-model/2735/2)) as follows:

1. Load the nrrd-file in slicer
  - This creates the `image` volume as expected

2. Create a model from the `image` volume
  - 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`

3. 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.

4. 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

Thanks for your help

---

<div class="post-metadata">

**Author:** ![lassoan](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/lassoan/32/13_2.png) [@lassoan](https://discourse.slicer.org/u/lassoan)\
**Post date:** [April 20, 2020, 2:32pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/2 "2020-04-20T14:32:51Z")

</div>

> [@hcmidt](#):
>
> Create a model from the `image` volume

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.

> [@hcmidt](#):
>
> The `probed_model` appears as a green sphere. All the surface is uniformly covered with green color.

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.

---

<div class="post-metadata">

**Author:** ![pieper](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/pieper/32/8_2.png) [@pieper](https://discourse.slicer.org/u/pieper)\
**Post date:** [April 20, 2020, 2:50pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/3 "2020-04-20T14:50:07Z")

</div>

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.

> <https://gist.github.com/pieper/f9da3e0a73c70981b48d0747132526d5>

---

<div class="post-metadata">

**Author:** ![hcmidt](https://avatars.discourse-cdn.com/v4/letter/h/59ef9b/32.png) [@hcmidt](https://discourse.slicer.org/u/hcmidt)\
**Post date:** [April 20, 2020, 5:51pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/4 "2020-04-20T17:51:54Z")

</div>

> [@lassoan](#):
>
> 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.

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.

> [@lassoan](#):
>
> 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 is a link with the dataset that results on my machine, following the steps I described above: [probed\_model.vtk - Google Drive](https://drive.google.com/open?id=1MN9jkQvmC0Ihu_MWEvcrCunLGFjcGTPe)

---

<div class="post-metadata">

**Author:** ![lassoan](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/lassoan/32/13_2.png) [@lassoan](https://discourse.slicer.org/u/lassoan)\
**Post date:** [April 20, 2020, 6:23pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/5 "2020-04-20T18:23:40Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![hcmidt](https://avatars.discourse-cdn.com/v4/letter/h/59ef9b/32.png) [@hcmidt](https://discourse.slicer.org/u/hcmidt)\
**Post date:** [April 20, 2020, 6:24pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/6 "2020-04-20T18:24:40Z")

</div>

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.

Running

```auto
modelNode = getNode('probed_model')
modelPointValues = modelNode.GetPolyData().GetPointData().GetArray("Normals")

```

in the slicer python prompt returns `None`. Running the same with the `sperical_model` returns a vtkFloatArray object that I can inspect.

So I suppose something might go wrong during the probing stage?

Thanks!

---

<div class="post-metadata">

**Author:** ![hcmidt](https://avatars.discourse-cdn.com/v4/letter/h/59ef9b/32.png) [@hcmidt](https://discourse.slicer.org/u/hcmidt)\
**Post date:** [April 20, 2020, 6:27pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/7 "2020-04-20T18:27:32Z")

</div>

> [@lassoan](#):
>
> 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?

Yes, here it is: [VolumeProperty.vp - Google Drive](https://drive.google.com/open?id=1AcxCfIc9wHuoVL2Q4yWQyDINp0LZR4o3)

---

<div class="post-metadata">

**Author:** ![lassoan](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/lassoan/32/13_2.png) [@lassoan](https://discourse.slicer.org/u/lassoan)\
**Post date:** [April 20, 2020, 6:28pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/8 "2020-04-20T18:28:57Z")

</div>

> [@hcmidt](#):
>
> in the slicer python prompt returns `None` . Running the same with the `sperical_model` returns a vtkFloatArray object that I can inspect.
> 
> So I suppose something might go wrong during the probing stage?

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).

---

<div class="post-metadata">

**Author:** ![lassoan](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/lassoan/32/13_2.png) [@lassoan](https://discourse.slicer.org/u/lassoan)\
**Post date:** [April 20, 2020, 6:30pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/9 "2020-04-20T18:30:10Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![hcmidt](https://avatars.discourse-cdn.com/v4/letter/h/59ef9b/32.png) [@hcmidt](https://discourse.slicer.org/u/hcmidt)\
**Post date:** [April 20, 2020, 6:40pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/10 "2020-04-20T18:40:50Z")

</div>

Perfect, here is the whole package: [https://drive.google.com/open?id=1OsToo4BLIQldFbSPRsOtUhPuqoM1vpX3](https://drive.google.com/open?id=1OsToo4BLIQldFbSPRsOtUhPuqoM1vpX3)

---

<div class="post-metadata">

**Author:** ![lassoan](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/lassoan/32/13_2.png) [@lassoan](https://discourse.slicer.org/u/lassoan)\
**Post date:** [April 20, 2020, 6:59pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/11 "2020-04-20T18:59:53Z")

</div>

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](https://www.slicer.org/wiki/Documentation/Nightly/ScriptRepository#Show_a_simple_surface_mesh_as_a_model_node) how to create a surface mesh using a VTK source. You can use a `vtkSphereSource` to create a sphere model.

---

<div class="post-metadata">

**Author:** ![hcmidt](https://avatars.discourse-cdn.com/v4/letter/h/59ef9b/32.png) [@hcmidt](https://discourse.slicer.org/u/hcmidt)\
**Post date:** [April 20, 2020, 7:37pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/12 "2020-04-20T19:37:25Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![lassoan](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/lassoan/32/13_2.png) [@lassoan](https://discourse.slicer.org/u/lassoan)\
**Post date:** [April 21, 2020, 12:17am UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/13 "2020-04-21T00:17:14Z")

</div>

> [@hcmidt](#):
>
> `Probe volume with model` module, the surface values would change according to the values of the volume

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.

> [@hcmidt](#):
>
> 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 did it well. Just don’t create your model by iso-surfacing the same volume that you use for coloring.

> [@hcmidt](#):
>
> Then back to the original problem: how can I achieve what I want to do in slicer?

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](https://vtk.org/doc/nightly/html/classvtkImageToPoints.html), remove zero-valued points using [thresholding](https://vtk.org/doc/nightly/html/classvtkThreshold.html), then create a surface using one of the [surface reconstruction filters](https://lorensen.github.io/VTKExamples/site/Cxx/Points/PoissonExtractSurface/). These all should not be more than 10-20 lines of Python code in total.

---

<div class="post-metadata">

**Author:** ![hcmidt](https://avatars.discourse-cdn.com/v4/letter/h/59ef9b/32.png) [@hcmidt](https://discourse.slicer.org/u/hcmidt)\
**Post date:** [April 21, 2020, 11:13am UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/14 "2020-04-21T11:13:04Z")

</div>

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.

Thank you!

---

<div class="post-metadata">

**Author:** ![hcmidt](https://avatars.discourse-cdn.com/v4/letter/h/59ef9b/32.png) [@hcmidt](https://discourse.slicer.org/u/hcmidt)\
**Post date:** [April 24, 2020, 2:07pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/15 "2020-04-24T14:07:46Z")

</div>

Hello together,

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):

```auto
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:

> **[example.nrrd](https://drive.google.com/file/d/1MoxMcTeA0RqUW67VGNldHtLc1qyZmTRL/view?usp=drive_open)**
>
> Google Drive file.

The workflow would be the following:

1. load `example.nrrd` in slicer
2. open the segment editor and extract the segment by applying the threshold operation.
3. run the `get_colored_model` function from the python interactor within slicer
4. 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!

---

<div class="post-metadata">

**Author:** ![hcmidt](https://avatars.discourse-cdn.com/v4/letter/h/59ef9b/32.png) [@hcmidt](https://discourse.slicer.org/u/hcmidt)\
**Post date:** [April 24, 2020, 8:07pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/16 "2020-04-24T20:07:12Z")

</div>

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’ll try to work out (a) next.

---

<div class="post-metadata">

**Author:** ![lassoan](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/lassoan/32/13_2.png) [@lassoan](https://discourse.slicer.org/u/lassoan)\
**Post date:** [April 24, 2020, 8:26pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/17 "2020-04-24T20:26:57Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![hcmidt](https://avatars.discourse-cdn.com/v4/letter/h/59ef9b/32.png) [@hcmidt](https://discourse.slicer.org/u/hcmidt)\
**Post date:** [April 24, 2020, 9:54pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/18 "2020-04-24T21:54:12Z")

</div>

Thanks! I was looking for this functionality before but did not manage to find it.  
Could you tell me how I can do that on the python level?

---

<div class="post-metadata">

**Author:** ![lassoan](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/lassoan/32/13_2.png) [@lassoan](https://discourse.slicer.org/u/lassoan)\
**Post date:** [April 24, 2020, 10:48pm UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/19 "2020-04-24T22:48:57Z")

</div>

You can find Python examples for all commonly needed operations and complete workflow examples in the [script repository](https://www.slicer.org/wiki/Documentation/Nightly/ScriptRepository#Segmentations).

---

<div class="post-metadata">

**Author:** ![hcmidt](https://avatars.discourse-cdn.com/v4/letter/h/59ef9b/32.png) [@hcmidt](https://discourse.slicer.org/u/hcmidt)\
**Post date:** [April 25, 2020, 2:14am UTC](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202/20 "2020-04-25T02:14:07Z")

</div>

Thanks. I have been looking at the repository and I am managing to extract the modelNode like so:

```auto
modelHierarchyNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLModelHierarchyNode')
slicer.modules.segmentations.logic().ExportVisibleSegmentsToModelHierarchy(segmentationNode, modelHierarchyNode)
modelNode = slicer.mrmlScene.GetFirstNodeByClass("vtkMRMLModelNode")

```

However, I _don’t manage to project the volume values onto the scalar values of the model like I did in my example_.

In case this is informative: the coordinates returned by  
`vtk_to_numpy(modelNode.GetPolyData().GetPoints().GetData())` are only four points

```auto
[[54.22350311, -200.5 , 100.],
 [-254.22349548, -200.5 , 100.],
 [54.22350311, 0.5 , 100.],
 [-254.22349548, 0.5 , 100.]]

```

which I found a bit strange.

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 :\>

[Next page](https://discourse.slicer.org/t/coloring-surfaces-in-slicer/11202.md?page=2)