Vedo: Setting pointColors with custom colormaps

[unidesigner]

The results of the mesh_coloring.py example does not produce a colored man in the middle as would be expected, see image:

What I'd like to do is setting the colors per point of an polydata object based on my own colormap, so that I can define the scaling from min to max.

# myarr is a numpy float32 1D array with e.g. min=10, max = 20000
from vtkplotter import colorMap
r, g, b = colorMap(... name='jet', vmin=myarr.min(), vmax=myarr.max())

and somehow pass this to an actor:

man1.pointColors(myarr alpha=0.5, cmap=MyColormapAbove)

A working example would be very helpful.

[marcomusy]

this is something that already happened to me on a laptop, I hope I can reproduce it, on my desktop I don't see the issue:

what vtk/python version do you have? OS?

[unidesigner]

ipython 3.6.7, vtk-8.1.2. on an ubuntu 18.04. I could try the vtk version shipped by ubuntu, or try anaconda just to see if it also happens there.

[marcomusy]

Your specs look fine..
Try out this, which is probably what you are aiming for:

from vtkplotter import Plotter

# these are the available color maps
mapkeys = ['afmhot', 'binary', 'bone', 'cool', 'coolwarm', 'copper', 
           'gist_earth', 'gray', 'hot', 'jet', 'rainbow', 'winter']

vp = Plotter(N=len(mapkeys), axes=4)

mug = vp.load('data/shapes/mug.ply')
scalars = mug.coordinates()[:,1] # let y-coord be the scalar

for i,key in enumerate(mapkeys): # for each available color map name
    imug = mug.clone()
    imug.pointColors(scalars, cmap=key)    
    vp.show(imug, at=i, legend=key)

vp.show(interactive=1)

if it doesn't work I need to see if I can reproduce the issue on my old laptop :)

[unidesigner]

thanks for the example code. unfortunately, i cannot make it work. i now also tested the code using miniconda3, doing a conda install vtk before doing a 'python3 setup.py install' to install from source. if you have other suggestions to try, i could look into it. or you try to reproduce it on your old laptop :)

[marcomusy]

It works fine also on my laptop now... ;(

what is the output of this for you?

from __future__ import print_function
import vtk

def MakeLUT(tableSize):
    '''
    Make a lookup table from a set of named colors.
    :param: tableSize - The table size
    :return: The lookup table.
    '''
    nc = vtk.vtkNamedColors()

    lut = vtk.vtkLookupTable()
    lut.SetNumberOfTableValues(tableSize)
    lut.Build()

    # Fill in a few known colors, the rest will be generated if needed
    lut.SetTableValue(0,nc.GetColor4d("Black"))
    lut.SetTableValue(1,nc.GetColor4d("Banana"))
    lut.SetTableValue(2,nc.GetColor4d("Tomato"))
    lut.SetTableValue(3,nc.GetColor4d("Wheat"))
    lut.SetTableValue(4,nc.GetColor4d("Lavender"))
    lut.SetTableValue(5,nc.GetColor4d("Flesh"))
    lut.SetTableValue(6,nc.GetColor4d("Raspberry"))
    lut.SetTableValue(7,nc.GetColor4d("Salmon"))
    lut.SetTableValue(8,nc.GetColor4d("Mint"))
    lut.SetTableValue(9,nc.GetColor4d("Peacock"))

    return lut

def MakeLUTFromCTF(tableSize):
    '''
    Use a color transfer Function to generate the colors in the lookup table.
    See: http://www.vtk.org/doc/nightly/html/classvtkColorTransferFunction.html
    :param: tableSize - The table size
    :return: The lookup table.
    '''
    ctf = vtk.vtkColorTransferFunction()
    ctf.SetColorSpaceToDiverging()
    # Green to tan.
    ctf.AddRGBPoint(0.0, 0.085, 0.532, 0.201)
    ctf.AddRGBPoint(0.5, 0.865, 0.865, 0.865)
    ctf.AddRGBPoint(1.0, 0.677, 0.492, 0.093)

    lut = vtk.vtkLookupTable()
    lut.SetNumberOfTableValues(tableSize)
    lut.Build()

    for i in range(0,tableSize):
        rgb = list(ctf.GetColor(float(i)/tableSize))+[1]
        lut.SetTableValue(i,rgb)

    return lut

def MakeCellData(tableSize, lut, colors):
    '''
    Create the cell data using the colors from the lookup table.
    :param: tableSize - The table size
    :param: lut - The lookup table.
    :param: colors - A reference to a vtkUnsignedCharArray().
    '''
    for i in range(1,tableSize):
        rgb = [0.0, 0.0, 0.0]
        lut.GetColor(float(i) / (tableSize - 1),rgb)
        ucrgb = list(map(int, [x * 255 for x in rgb]))
        colors.InsertNextTuple3(ucrgb[0], ucrgb[1], ucrgb[2])
        s = '['+ ', '.join(['{:0.6f}'.format(x) for x in rgb]) + ']'
        print(s, ucrgb)


def main():
    '''
    :return: The render window interactor.
    '''

    nc = vtk.vtkNamedColors()

    # Provide some geometry
    resolution = 3

    plane11 = vtk.vtkPlaneSource()
    plane11.SetXResolution(resolution)
    plane11.SetYResolution(resolution)

    plane12 = vtk.vtkPlaneSource()
    plane12.SetXResolution(resolution)
    plane12.SetYResolution(resolution)

    tableSize = max(resolution * resolution + 1, 10)

    #  Force an update so we can set cell data
    plane11.Update()
    plane12.Update()

    #  Get the lookup tables mapping cell data to colors
    lut1 = MakeLUT(tableSize)
    lut2 = MakeLUTFromCTF(tableSize)

    colorData1 = vtk.vtkUnsignedCharArray()
    colorData1.SetName('colors') # Any name will work here.
    colorData1.SetNumberOfComponents(3)
    print('Using a lookup table from a set of named colors.')
    MakeCellData(tableSize, lut1, colorData1)
    # Then use SetScalars() to add it to the vtkPolyData structure,
    # this will then be interpreted as a color table.
    plane11.GetOutput().GetCellData().SetScalars(colorData1)

    colorData2 = vtk.vtkUnsignedCharArray()
    colorData2.SetName('colors') # Any name will work here.
    colorData2.SetNumberOfComponents(3)
    print('Using a lookup table created from a color transfer function.')
    MakeCellData(tableSize, lut2, colorData2)
    plane12.GetOutput().GetCellData().SetScalars(colorData2)

    # Set up actor and mapper
    mapper11 = vtk.vtkPolyDataMapper()
    mapper11.SetInputConnection(plane11.GetOutputPort())
    # Now, instead of doing this:
    # mapper11.SetScalarRange(0, tableSize - 1)
    # mapper11.SetLookupTable(lut1)
    # We can just use the color data that we created from the lookup table and
    # assigned to the cells:
    mapper11.SetScalarModeToUseCellData()
    mapper11.Update()

    mapper12 = vtk.vtkPolyDataMapper()
    mapper12.SetInputConnection(plane12.GetOutputPort())
    mapper12.SetScalarModeToUseCellData()
    mapper12.Update()

    writer = vtk.vtkXMLPolyDataWriter()
    writer.SetFileName('pdlut.vtp')
    writer.SetInputData(mapper11.GetInput())
    # This is set so we can see the data in a text editor.
    writer.SetDataModeToAscii()
    writer.Write()
    writer.SetFileName('pdctf.vtp')
    writer.SetInputData(mapper12.GetInput())
    writer.Write()

    actor11 = vtk.vtkActor()
    actor11.SetMapper(mapper11)
    actor12 = vtk.vtkActor()
    actor12.SetMapper(mapper12)

    # Let's read in the data we wrote out.
    reader1 = vtk.vtkXMLPolyDataReader()
    reader1.SetFileName("pdlut.vtp")

    reader2 = vtk.vtkXMLPolyDataReader()
    reader2.SetFileName("pdctf.vtp")

    mapper21 = vtk.vtkPolyDataMapper()
    mapper21.SetInputConnection(reader1.GetOutputPort())
    mapper21.SetScalarModeToUseCellData()
    mapper21.Update()
    actor21 = vtk.vtkActor()
    actor21.SetMapper(mapper11)

    mapper22 = vtk.vtkPolyDataMapper()
    mapper22.SetInputConnection(reader2.GetOutputPort())
    mapper22.SetScalarModeToUseCellData()
    mapper22.Update()
    actor22 = vtk.vtkActor()
    actor22.SetMapper(mapper22)

    # Define viewport ranges.
    # (xmin, ymin, xmax, ymax)
    viewport11 = [0.0, 0.0, 0.5, 0.5]
    viewport12 = [0.0, 0.5, 0.5, 1.0]
    viewport21 = [0.5, 0.0, 1.0, 0.5]
    viewport22 = [0.5, 0.5, 1.0, 1.0]

    # Set up the renderers.
    ren11 = vtk.vtkRenderer()
    ren12 = vtk.vtkRenderer()
    ren21 = vtk.vtkRenderer()
    ren22 = vtk.vtkRenderer()

    # Setup the render windows
    renWin = vtk.vtkRenderWindow()
    renWin.SetSize(800, 800)
    renWin.AddRenderer(ren11)
    renWin.AddRenderer(ren12)
    renWin.AddRenderer(ren21)
    renWin.AddRenderer(ren22)
    ren11.SetViewport(viewport11)
    ren12.SetViewport(viewport12)
    ren21.SetViewport(viewport21)
    ren22.SetViewport(viewport22)
    ren11.SetBackground(nc.GetColor3d('MidnightBlue'))
    ren12.SetBackground(nc.GetColor3d('MidnightBlue'))
    ren21.SetBackground(nc.GetColor3d('MidnightBlue'))
    ren22.SetBackground(nc.GetColor3d('MidnightBlue'))
    ren11.AddActor(actor11)
    ren12.AddActor(actor12)
    ren21.AddActor(actor21)
    ren22.AddActor(actor22)

    iren = vtk.vtkRenderWindowInteractor()
    iren.SetRenderWindow(renWin)
    renWin.Render()

    return iren

def CheckVTKVersion(requiredMajorVersion):
    '''
    Check the VTK version.
    :param: requiredMajorVersion e.g. 6
    '''
    version = vtk.vtkVersion()
    if version.GetVTKMajorVersion() > requiredMajorVersion:
        raise
    else:
        return

if __name__ == '__main__':
    try:
        pass
    except:
        print("You need VTK Version 6 or greater.")
        print("The class vtkNamedColors is in VTK version 6 or greater.")
        exit(0)
    iren = main()
    iren.Start()

[unidesigner]

the output is the same as you have

[marcomusy]

..what about this:

import vtk
from vtkplotter import Plotter, colors
import numpy as np
from vtk.util.numpy_support import numpy_to_vtk

vp = Plotter(verbose=0)

g = vp.grid().wire(False)

scalars = g.coordinates()[:,0] # let x-coord be the scalar

# g.pointColors(scalars, cmap='jet')

vmin, vmax = np.min(scalars), np.max(scalars)
lut = vtk.vtkLookupTable()
lut.SetTableRange(vmin, vmax)
n = len(scalars)
lut.SetNumberOfTableValues(n)
lut.Build()

for i in range(n):
    c = colors.colorMap(i, 'jet', 0, n)
    lut.SetTableValue(i, c[0], c[1], c[2], 1)

arr = numpy_to_vtk(np.ascontiguousarray(scalars), deep=True)
arr.SetName('pointColors_jet')
print([arr.GetNumberOfComponents()])
mapper = g.GetMapper()
mapper.SetScalarRange(vmin, vmax)
mapper.SetLookupTable(lut)
mapper.ScalarVisibilityOn()
vtkpts = g.polydata().GetPointData()
# vtkpts.AddArray(arr)
vtkpts.SetScalars(arr)
vtkpts.SetActiveScalars('pointColors_jet')

vp.show()

(thanks for your patience for debugging this..)

[unidesigner]

It's an interesting visual git bisect exercise :)
In this case, the square is grey for me.

[marcomusy]

uhmm, I suspect of numpy_to_vtk .. what about this:

import vtk
from vtkplotter import Plotter, colors
import numpy as np
from vtk.util.numpy_support import numpy_to_vtk

vp = Plotter(verbose=0)

g = vp.grid().wire(False)

scalars = g.coordinates()[:,0] # let x-coord be the scalar

# g.pointColors(scalars, cmap='jet')

vmin, vmax = np.min(scalars), np.max(scalars)
lut = vtk.vtkLookupTable()
lut.SetTableRange(vmin, vmax)
n = len(scalars)
lut.SetNumberOfTableValues(n)
lut.Build()

for i in range(n):
    R,G,B = colors.colorMap(i, 'jet', 0, n)
    lut.SetTableValue(i, R,G,B, 1)

dummyarr = numpy_to_vtk(np.ascontiguousarray(scalars), deep=True)
print('numpy_to_vtk, nr. of tuples', dummyarr.GetNumberOfTuples())

arr = vtk.vtkFloatArray()
for s in scalars:
    arr.InsertNextTuple1(s)
arr.SetName('pointColors_jet')
print('vtkFloatArray, nr. of tuples', arr.GetNumberOfTuples())

mapper = g.GetMapper()
mapper.SetScalarRange(vmin, vmax)
mapper.SetLookupTable(lut)
mapper.ScalarVisibilityOn()
vtkpts = g.polydata().GetPointData()
# vtkpts.AddArray(arr)
vtkpts.SetScalars(arr)
vtkpts.SetActiveScalars('pointColors_jet')

vp.show()

what is your print output?

PS: are you sure you can import matplotlib.cm

[unidesigner]

thanks for your patient. as it turns out, the problem was indeed that I could not import matplotlib.cm
doing a sudo apt-get install python3-matplotlib solved all the issues and i can see the proper colormappings in all of the examples above! great that this works now.

perhaps an exception should be thrown somewhere when using colormaps related function calls and matplotlib is not installed.

[marcomusy]

Fantastic! I've added a exception message for matplotlib.
In the process I also found another bug which shows up when the minimum scalar value is not zero.
Both will be fixed in the next commit.

[marcomusy]

I think I fixed a couple of other issues.. and committed a new version. I also added an other similar example here: https://github.com/marcomusy/vtkplotter/blob/master/examples/basic/mesh_bands.py

[unidesigner]

Thanks, this is very useful to see how to adjust alpha values with an array.

One follow-up question: I add multiple mesh actors to a scene with their respective per-point scalar values. Now I'd like to only show one colormap that applies to all meshes equally, irrespective of their own range of per-point scalar values. How would this be possible with the current API and not diving into VTK myself deeper? E.g. I do not quite understand how vtkplotter.colorMap (which would have vmin/vmax) can be applied to individual meshes. Thanks for help.

[marcomusy]

Hi, you can just use for all meshes the same vmin and vmax, like in this example:

in the first renderer the scalar runs from 28 to 44, in the second from 18 to 34, but they share a unique color mapping.

from vtkplotter import Plotter

vp = Plotter(N=2)

##################################### 
man1 = vp.load('data/shapes/man.vtk')
scals = man1.coordinates()[:,2]*5 + 37 # pick z coordinates of vertices [28->44]

man1.pointColors(scals, cmap='jet', vmin=18, vmax=44) 
vp.show(man1, at=0, axes=0)

##################################### 
man2 = vp.load('data/shapes/man.vtk')
scals = man2.coordinates()[:,2]*5 + 27 # pick z coordinates of vertices [18->34]

man2.pointColors(scals, cmap='jet', vmin=18, vmax=44) 
vp.show(man2, at=1, axes=0)

vp.addScalarBar()
vp.show(interactive=1)

(if the output shows different I will do a commit of my local version)

[unidesigner]

I get a TypeError: pointColors() got an unexpected keyword argument 'vmin'

[marcomusy]

uhmm, are you sure you're using the latest version (8.8.0) ?

[unidesigner]

Indeed, there was an old version hiding which got imported first. Your example works now and I can control the colormap with the vmin/vmax as I wanted. Thanks.

[NaGenhao]

How to set edgeColors by it's edge length?

L = Lines(startPoints, endPoints)
L.cellColors(len_edges)

[marcomusy]

sorry that's not very clear to me, do you want to set colors to cell based on edge length, or color the edges?
Is the snippet already the solution? (Indeed it is a very clever one to colorize lines!)

from vtkplotter import *
import numpy as np

sph = Sphere().cutWithPlane()

sPoints = np.zeros([sph.N(),3]) + [0.75, 0, 0.1]
ePoints = sph.coordinates()

rays = Lines(sPoints, ePoints)
rays.cellColors(mag(sPoints-ePoints), cmap="spring")
rays.show(bg="w", axes=2)