# How to colorize surface by slope?

I’ve subclassed vtkAbstractPolyDataReader to MyReader, which reads topographic data from a file and outputs vtkPolyData containing points and polygons (triangles). I’ve written an app that uses MyReader to read a topographic data file and display it as a 3D surface, colorized by elevation, like so:

// Read topographic data

// Want to colorize topographic data by elevation
vtkSmartPointer<vtkElevationFilter> colorizer =
vtkSmartPointer<vtkElevationFilter>::New();



[and so on]

For example, the first image (below) is colored by elevation (depth) of Monterey Canyon - blue is deep, yellow-orange are shallow.

Now I want my app to colorize by slope gradients, as shown in the second image of Monterey Canyon; red are the steepest slopes, blue (on the canyon floor) is nearly flat.

I think I can do this somehow with vtkGradientFilter, but I’m not sure where to begin. The few examples I’ve seen - VTKExamples Gradient and GradientFilter are not very clear to me. Could anyone advise on how I would change my app to achieve this?
Thanks!

Hello,

Maybe you could post your entire scene-building code so one can easily identify where in you current pipeline to insert a vtkGradientFilter.

regards,

Paulo

1 Like

I’ve recently created an example for vedo:
vedo -r isolines

relevant lines of python code are:

        on='cells'
if on.startswith('p'):
varr = self.inputdata().GetPointData()
tp = vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS
else:
varr = self.inputdata().GetCellData()
tp = vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS
gra.SetInputData(inputdata)
gra.SetInputScalars(tp, arrname)
gra.Update()
return gra.GetOutput()

1 Like

Hi @Paulo_Carvalho - here is a more complete listing of my code, that currently colors by elevation with vtkElevationFilter:

   // Read the grid file

std::cout << "*** reader->Update()" << std::endl;

std::string filePath = argv[argc-1];
if (filePath[0] == '-') {
std::cerr << filePath << ": invalid grid file name" << std::endl;
return EXIT_FAILURE;
}

std::cout << "*** reader->Update()" << std::endl;

// Color data points based on z-value
vtkSmartPointer<vtkElevationFilter> elevation =
vtkSmartPointer<vtkElevationFilter>::New();

float zMin, zMax;
elevation->SetLowPoint(0, 0, zMin);
elevation->SetHighPoint(0, 0, zMax);

// Visualize the data...

// Create renderer
vtkSmartPointer<vtkRenderer> renderer =
vtkSmartPointer<vtkRenderer>::New();

// Create mapper
vtkSmartPointer<vtkPolyDataMapper> mapper =
vtkSmartPointer<vtkPolyDataMapper>::New();

mapper->SetInputConnection(elevation->GetOutputPort());

if (useLUT) {
// Use lookup table
elevation->SetScalarRange(zMin, zMax);

// Create the lookup table
vtkSmartPointer<vtkLookupTable> lut = vtkSmartPointer<vtkLookupTable>::New();
// Fill the lookup table
makeLUT(colorScheme, lut);

mapper->SetScalarRange(zMin, zMax);
mapper->ScalarVisibilityOn();
mapper->SetLookupTable(lut);
}

// Create actor
vtkSmartPointer<vtkActor> actor =
vtkSmartPointer<vtkActor>::New();

// Assign mapper to actor
actor->SetMapper(mapper);

// Add actor to the renderer

// Create renderWindow
vtkSmartPointer<vtkRenderWindow> renderWindow =
vtkSmartPointer<vtkRenderWindow>::New();

// Add renderer to the renderWindow

// Create renderWindowInteractor
vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
vtkSmartPointer<vtkRenderWindowInteractor>::New();

vtkSmartPointer<vtkInteractorStyleTrackballCamera> style =
vtkSmartPointer<vtkInteractorStyleTrackballCamera>::New();

renderWindowInteractor->SetInteractorStyle(style);

renderWindowInteractor->SetRenderWindow(renderWindow);

Hi, Tomasso,

I’ve made changes to your pipeline (signaled with a //<-------------------). I didn’t have time to test it, so read it to see whether it makes sense for you, then give it a try. It’s likely it won’t work in your first attempt, but that’s a kickstart.

regards,

Paulo

1 Like

Hi @Paulo_Carvalho - thanks very much for your suggestions! I incorporated them into my code. When I run my app I see about a 15 second delay before the surface is displayed, presumably because the gradient filter is computing slopes (dataset consists of 1113 x 1450 points). However the final displayed surface colors are identical to elevation colors, i.e. they don’t reflect the gradient. Why might that be?
Also, if I want to color the surface by gradient and not by elevation, why does my pipeline need to include an elevation filter?

Here is the current code:

  vtkSmartPointer<mb_system::GmtGridReader> reader =

float zMin, zMax;

// Color data points based on z-value
vtkSmartPointer<vtkElevationFilter> elevationFilter =
vtkSmartPointer<vtkElevationFilter>::New();

elevationFilter->SetLowPoint(0, 0, zMin);
elevationFilter->SetHighPoint(0, 0, zMax);

// Visualize the data...
// Create renderer
vtkSmartPointer<vtkRenderer> renderer =
vtkSmartPointer<vtkRenderer>::New();

// Create mapper
vtkSmartPointer<vtkPolyDataMapper> mapper =
vtkSmartPointer<vtkPolyDataMapper>::New();

if (useLUT) {

vtkSmartPointer<vtkLookupTable> lut =
vtkSmartPointer<vtkLookupTable>::New();

makeLUT(colorScheme, lut);

mapper->SetScalarRange(zMin, zMax);
mapper->ScalarVisibilityOn();
mapper->SetLookupTable(lut);
}

// Create actor
vtkSmartPointer<vtkActor> actor =
vtkSmartPointer<vtkActor>::New();

// Assign mapper to actor
actor->SetMapper(mapper);

// Create renderWindow
vtkSmartPointer<vtkRenderWindow> renderWindow =
vtkSmartPointer<vtkRenderWindow>::New();

// Add renderer to the renderWindow
renderWindow->AddRenderer(renderer);

Hi, Tom,

Attack one problem at a time. It seems that the vtkDataSet in the output of the gradient filter has two or more scalars and the first field (the default one) is probably the one with elevation values. You have to change your pipeline to use the one with the gradient values (e.g. with gradientFilter->GetOutput()-> [ GetPointData() | GetCellData() ] ->SetActiveScalars('myScalarField').

regards,

Paulo

Thank you @Paulo_Carvalho - sorry , I really am a novice and I really appreciate your advice. How do I know what data is available from the gradientFilter, i.e. how many scalar fields, and how do I know to choose gradientFilter’s GetPointData() versus GetCellData()? When I’ve made that choice then from other examples it seems the pattern is this:

vtkSmartPointer<vtkGradientFilter> gradientFilter =

// Make sure mapper uses correct data

[ . . .]

vtkSmartPointer<vtkPolyDataMapper> mapper =
vtkSmartPointer<vtkPolyDataMapper>::New();

mapper->SetInputConnection(gradientFilter->GetOutputPort());

vtkGradientFilter's GetOutput() returns a pointer to a vtkDataSet. A vtkDataSet is a very generic class than can be a point cloud, a polygonal line, a grid, a tetrahedron mesh, you name it. Here’s the documentation of vtkDataSet: https://vtk.org/doc/nightly/html/classvtkDataSet.html . See how many subclasses it has. Quite a lot! We can assume that the filter returns a data set with the same sub-type as the input’s. This is just to introduce you to the basics of VTK design. It’s not really relevant to your problem.

Now to the difference between GetPointData() versus GetCellData(). The former is a pointer to the data values stored in the points of your object and the latter points to the values stored in “cells”. The points are like “atoms” in XYZ space. A “cell” is a construct of “atoms”. A “cell” is made of one or more “atoms”. “Cells” can share “atoms” to optimize de model. A vtkVertex is a “cell” made with just one point. So you can say a point cloud is made of cells of type vtkVertex, each defined by a single point. Other types of “cells” are lines, triangles, cubes, etc. So, the data can reside directly in the points and/or in the cells. You have to know your input data set or simply test it to see where your reader is putting the data values.

For how many fields your data set has, you can try something like this:

vtkDataSet* dataSet = gradientFilter->GetOutput();
vtkCellData* cellData = dataSet->GetCellData();
vtkDataArray* dataArray = cellData->GetScalars();
std::cout << "My data set's cells have " << dataArray->GetNumberOfComponents() << " scalar field(s)." << std::endl;


I hope this helps to cast some light on VTK data sets.

cheers,

Paulo

2 Likes

just

pip install vedo
vedo -r isolines

1 Like

Thanks very much @marcomusy, this code is very useful. I’ve implemented your python gradient() function in C++, but I don’t quite understand how this section of isolines.py works - in particular, how exactly do you put the gradient colormap directly on the hillside mesh? Can you describe that process strictly in terms of the VTK C++ library? Many thanks!

gvecs = mesh.gradient('Scalars', on='cells')
cc = mesh.cellCenters()
ars = Arrows(cc, cc + gvecs*0.01, c='bone_r').lighting('off')
show(mesh, isol, ars, "Arrows=\nablaH", at=2)

# colormap the gradient magnitude directly on the mesh
mesh2 = mesh.clone().lw(0.1).cmap('jet', mag(gvecs), on='cells')
cam = dict(pos=(2.57, -1.92, 3.25), # get these nrs by pressing C
focalPoint=(0.195, 1.06, 8.65e-3),
viewup=(-0.329, 0.563, 0.758),
distance=5.00,
clippingRange=(1.87, 8.79))
show(mesh2, "Color=|\nablaH|", at=3, camera=cam, interactive=True)

Once you have the vtk array, you create a vtkLookupTable which defines the color associated to the scalar and pass it to the mapper; pass the array to either cells or vertices (see detailed explanation by Paulo). Something like

        lut = vtk.vtkLookupTable()
lut.SetNumberOfTableValues(ncols)
for i, c in enumerate(...):
lut.SetTableValue(i, r, g, b, alpha)
lut.Build()

mapper.SetLookupTable(lut)
mapper.ScalarVisibilityOn()
mapper.SetScalarRange(vmin, vmax)
#poly.GetCellData().SetScalars(arr)
poly.GetCellData().SetActiveScalars(arrayName)
#poly.GetCellData().Modified()

2 Likes

Sorry, I’m still struggling a bit…
In your pseudocode, “mapper” holds the poly data for the x,y,z terrain points, right? Where is the terrain gradient data generated by vtkGradientFilter.GetOutput() in this code - is it “poly”? I am still missing how you overlay the gradient colors onto the 3D terrain surface.

Thanks!

mapper holds the vtkMapper and poly the vtkPolyData. Have look at this example too: https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/ColorCells/

So when you say poly is the vtkPolyData, you mean poly is the vtkGradientFilter output? I don’t see where poly is connected to the mapper in the pseudo code though (sorry, I’m sounding dumb even to myself at this point ). Really appreciate your help!

So when you say poly is the vtkPolyData, you mean poly is the vtkGradientFilter output?

Yes

I don’t see where poly is connected to the mapper in the pseudo code though

in the vtk example it’s all done here:

// Setup actor and mapper
vtkSmartPointer<vtkPolyDataMapper> mapper =
vtkSmartPointer<vtkPolyDataMapper>::New();
mapper->SetInputConnection(aPlane->GetOutputPort()); # <--
mapper->SetScalarRange(0, tableSize - 1);
mapper->SetLookupTable(lut);

1 Like