First steps in 3D visualizations that are worth a thousand words - a 30 min introduction to visualization.
The NeSI Team.
What you will learn:
- The difference between 2D and 3D visualization
- What makes a 3D visualization truly 3D
- actors
- light
- background
- How to visualize grids
- How to show data on grids
- Playing with colours
- Where to go from here
Visualization is a process for representing objects or information in an intuitive way, typically using our innate concepts for perspective, depth of field, colour, shading and lighting.
Both plotting and visualization aim to convey as much information as possible. The two approaches are complementary: plotting represents objects and information on a flat sheet of paper whereas visualization will show those objects in a 3D scene. In a visualization the viewer will typically interact with the object (through rotation of the camera, adjusting the light, etc.). Plotting tends to be static, hence it is more suitable for publishing.
For custom visualizations we recommend VTK, which is free and precompiled versions are easy to install. To run the examples we recommend to use to create a virtual environment
# Path to the virtual environment
VTKENV=...
python -m venv $VTKENV
(on some systems you may need to type python3 -m venv $VTKENV
,
activate the environment
source $VTKENV/bin/activate
and install VTK using
pip install vtk
To exit the environment, type deactivate
.
Many users may also want to consider VisIt or Paraview as these tools will require no or minimal programming. If your file format is support by any of these tools then you will be able to explore your data more quickly. Both these tools are free and binaries exist for Windows, Linux and Mac OS X.
Visualizations are made of elements that are assembled into a pipeline, essentially a workflow. Typically, a pipeline will involve reading data, applying filters, creating actors and rendering the actors in a scene. Filters extract information from the data; e.g. a contour filter might extract the surface for which a field satisfy a constraint. Actors are objects that be rendered; triangles, lines, points, etc. At the end of the day every actor is made of collections of triangles, lines and points.
At a minimum, a pipeline will have a source (reader, object, etc.), a mapper to draw and colour an actor.
source -> mapper -> actor
Visualization is a little like a theater play; there is a background, there are lights and actors. We'll start with a cone playing the role of an actor.
Try:
python scene/cone.py
(This may take some time as a lot of shared libraries get loaded.)
Next we'll refine the cone and add some lights, each with its own colour.
Try:
python scene/coneWithLight.py
Questions:
- How can I change the colour of the light?
- How can add another light?
More likely, your visualizations will involve gridded data so we'll need to understand how grids are represented in VTK. The main types of grid are structured and unstructructed. Structured grids have regular topology with arbitrary points. This means that given a set of indices describing a point we can always find the neighbours without additional information beyond the size of the grid. Unstructured grids build on top of points and cells and the arrangement can be arbitrary, i.e. one has to specify the topology. You can mix cells of different types in VTK, for instance hexahedra with prisms, lines and points. See [http://www.vtk.org/data-model/]for a full list of supported cells. Naturally, unstructured grids are more flexible than structured grids. The latter are however easier to use (and more efficient).
We'll start with a polar grid represented as a structured grid.
Try:
python grid/polar.py
Note how the structured grid is constructed. First, the set of coordinates is computed and stored in object coords
,
an array of doubles
coords.SetNumberOfComponents(3)
coords.SetNumberOfTuples(numPoints)
k = 0
for j in range(nrho):
rho = j * drho
for i in range(nthe):
the = i * dthe
x = rho * cos(the)
y = rho * sin(the)
z = 0.0
coords.SetTuple(k, (x, y, z))
k += 1
Then an object pts
of type vtkPoints is created with the points set to array coords
. Finally,
the structured grid object grid
sets the points (pts
).
Using a for loop in python can be a little slow for large grids. Here is a little known trick - it is possible to create all the data using the python numpy module and then pass the data directly to VTK:
thes = numpy.linspace(0., 2*numpy.pi, nthe)
rhos = numpy.linspace(0., 1., nrho)
tthes, rrhos = numpy.meshgrid(thes, rhos)
xyz = numpy.zeros((numPoints, 3), numpy.float64)
rrhos = numpy.reshape(rrhos, (numPoints,))
tthes = numpy.reshape(tthes, (numPoints,))
xyz[:, 0] = rrhos*numpy.cos(tthes)
xyz[:, 1] = rrhos*numpy.sin(tthes)
coords.SetNumberOfComponents(3)
coords.SetNumberOfTuples(numPoints)
coords.SetVoidArray(xyz, 3*numPoints, 1)
Finally it should be noted that all VTK objects work with 3D grids and coordinates. Because our grid is 2D we must set one of the dimensions to one:
grid.SetDimensions(nthe, nrho, 1)
Question:
- How can I change the resolution of the grid?
In VTK data don't generally exist without a grid but a grid can exist without data. So far we just built a grid but neglected to attach data to the grid. There are two types of data: point data and cell data. Point data belong the grid nodes whereas cell data associate to grid cells.
To set the point data use:
grid.GetPointData().SetScalars(data)
To set cell data use
grid.GetCellData().SetScalars(data)
We'll start with point data. Try:
python grid/polarWithPointData.py
Compare this to cell data. Try:
python grid/polarWithCellData.py
You'll notice that the cell get a solid colour in the case of cell data. Point data are interpolated linearly from each node to the neighbour node.
The default colour map in VTK maps low values to red and high values to blue! We need to fix this by adding a lookup table and a color bar actor.
A lookup table
lut = vtk.vtkLookupTable()
maps a data value to a colour, encoded as four floats in the range of 0 to 1 (red, green, blue and opacity). Here we'll set the colour range to vary from white to green, to light blue, red, purple and black.
ncolors = 64
lut.SetNumberOfColors(ncolors)
for i in range(ncolors):
x = 0.5*i*numpy.pi/(ncolors - 1.)
r = numpy.cos(3*x)**2
g = numpy.cos(1*x)**2
b = numpy.cos(5*x)**2
a = 1.0 # opacity
lut.SetTableValue(i, r, g, b, a)
lut.SetTableRange(fMin, fMax)
You will also need to tell the mapper to use the look up table
dataMapper.SetUseLookupTableScalarRange(1)
dataMapper.SetLookupTable(lut)
The color bar is an actor
cbar = vtk.vtkScalarBarActor()
cbar.SetLookupTable(lut)
which must be added to the scene
ren.AddActor(cbar)
Try:
python grid/polarWithPointDataLut.py
Question:
- What is the effect of changing the number of colours (ncolors)?
So far we worked in flat land -- when generating grid points we set the elevation z
to zero. Let's set the elevation of the value of another function:
rr2 = (xx**2 + yy**2).reshape((numPoints,))
xyz[:, 2] = 0.2*numpy.sin(10.*rr2)/numpy.sqrt(2.*rr2)
and run:
python grid/polarWithPointDataLutBump.py
The previous example was not really in 3D - it displayed a surface which is a 2D object.
Try:
python 3d/cube.py
to run a case with 3D data. We're facing a new problem: we only see the exterior of the cube.
One possibility is generate iso-surfaces from the data. Let's say we want to show the surfaces for which the data take values 0.2, 0.4 and 0.7:
contour = vtk.vtkContourFilter()
contour.SetNumberOfValues(3)
contour.SetValue(0, 0.2)
contour.SetValue(1, 0.4)
contour.SetValue(2, 0.7)
contour.SetInputData(grid)
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(contour.getOutputPort())
Try:
python 3d/cubeIsoSurface.py
Another possibility is to cut through the data. We'll create a vtkCutter
object and let VTK know that we will cut the domain with a vtkPlane
:
plane = vtk.vtkPlane()
plane.SetPosition(0.5, 0.5, 0.5)
plane.SetNormal(1., 0.2, -0.3)
knife = vtk.vtkCutter()
knife.SetCutFunction(plane)
Try:
python 3d/cubeCut.py
In the previous step we cut through our data and needed provide a plane with position and normal vectors to indicate where to cut. In most instance you will likely want to vary the cut plane. This can be done by hooking VTK up to a graphical user interface (for instance Qt). Another way to achieve this is by saving your data in a file
writer = vtk.vtkStructuredGridWriter()
writer.SetFileName('cube.vtk')
writer.SetInputData(grid)
writer.Update()
and use VisIt to explore your data.
Hope you enjoyed the tutorial! Feel free to fork and create pull-requests.