Voxel plotting with gnuplot 5.4

In this followup to our coverage of the release of gnuplot 5.4, we look more deeply at one of the new features: voxel plots. We only briefly touched on these plots in that article, but they are the most conspicuous addition in this release of the free-software graphing tool. Voxel plotting provides multiple ways to visualize 3D data, so it is worth looking at this new plot type in more detail.

Voxel grids

Previously, we introduced the simple physical system of an electric dipole, with two equal but opposite charges placed on the z axis at 0.25 and 0.75. It should be instructive to stick with this one problem, and explore different ways to visualize the 3D potential field produced by the charges.

Before diving in to voxel plotting proper, we need to understand that the familiar splot command has new 3D powers. splot used to stand for "surface plot", but can now do much more. For example, we can use it to plot the positions of our two charges like this:

    set xrange [0:1]; set yrange [0:1]; set zrange [0:1]
    set view 65,40
    set xyplane at -0.1
    set border 4095
    unset key
    $charges << EOD
    0 0 0.75 1
    0 0 0.25 -1
    EOD
    splot $charges using 1:2:3:4 with points\
        pointtype 7 pointsize 5 linecolor palette

This script produces this picture:

The first six lines of the script set the ranges of the display bounding box, the angle of view, the position of the bottom plane, and set the borders to surround the box on all sides. The next line, beginning with $charges, defines a "data block" consisting of the following two lines. Each line contains x, y, z, coordinates and, in the fourth column, the magnitude of the charge. The final command, broken over two lines, plots the two charges using their positions, extracted with the using 1:2:3 piece, and the charge value from the fourth column, extracted with the :4. This value is used to decide which colors the plotted points should be, by mapping the value onto the color palette, which is what the "linecolor palette" tells gnuplot to do. The other clauses set the pointsize to be five character widths and the pointtype to a circle (7).

Next, we will make a graph of the 3D structure of the potential field around these two charges. For this, we turn to the voxel grid. Just as a 2D image, such as a photograph, is a rectangular array of pixels, data in 3D can be represented as a 3D rectangular array of voxels, or volume pixels. Each voxel has x, y, and z coordinates, and a numerical value attached to it, so the voxel grid can represent a function of three variables, f(x, y, z). Note that this is completely new in gnuplot 5.4; previously, 3D plotting was confined to the plotting of surfaces or other representations of functions of two variables.

The way that the voxel commands in gnuplot are structured fits naturally with concepts from physics, which is why we've chosen a physical system for our example. In classical physics one has sources and the fields that they give rise to. The source could be a mass, such as a planet, creating a gravitational field; or, as in our case, a charged particle, giving rise to an electrical field. The field from a collection of sources is the sum of the fields from each individual source. The sources in our examples will be the two charges depicted above. The field created by adding the potentials from the two charges will be calculated and stored on the voxel grid.

Voxel plotting proceeds in two stages. First, after the voxel grid is initialized, it is filled with values calculated from any number of sources. This stage attaches the numbers to the locations on the voxel grid; it defines the f(x, y, z). In the second stage, we use the various styles of the splot command to create visualizations of this 3D function. Thus the calculation that establishes the f(x, y, z) need only be done once; the data is permanently stored in the voxel grid. Then we can make any number of plots of different types, from different angles, and so on. This two-stage process may seem cumbersome, but, since the initial filling of the voxel grid can be computationally intensive, it is actually more efficient to be able to do this just once.

There can be any number of sources. In our example, there are just two: the positive and negative charges. They need not correspond to voxel locations, and there can be more sources than voxel grid elements. The collection of sources, in other words, is independent of the configuration of the voxel grid.

The new command set vgrid establishes the voxel grid. It requires a name beginning with a dollar sign and a size; an example of a complete command would be:

    set vgrid $v size 25

which defines a 25 × 25 × 25 grid initialized with 0 at each voxel. Voxel grids only come in the cubic variety. This command also sets the voxel grid $v to be the currently active grid, so that certain subsequent commands will automatically operate on it.

The command to fill the voxel grid with values from the sources is vfill. A complete vfill command has this form:

    vfill source using x:y:z:radius:value

where source is a file or data block that must have at least three columns; in our example, that's the data block called $charges. Each line of source will become a source of data for the voxel grid. The first three colon-separated fields after the using keyword refer to the three spatial coordinates; if we put 1:2:3 here then they will be taken from the first three respective columns in source. The radius is required and indicates that the source values are to be applied on the voxel grid within a sphere with the indicated radius for each point in the source. The value indicates what to add to the voxel grid points inside that radius for each source point. This can be a constant number, some function of the coordinates, values taken from other optional columns in source, or anything else.

For example, say there was a file of 3D data, called "f", stored on disk; it has a 100 × 100 × 100 array, with x, y, z coordinates in the first three columns, ranging from 0 to 1 along each dimension, and the data value in the fourth column. This data could be assigned to a voxel grid directly with the following commands:

    set xrange [0 : 1]; set yrange [0 : 1]; set zrange [0 : 1]
    set vxrange [0 : 1]; set vyrange [0 : 1]; set vzrange [0 : 1]
    set vgrid $v size 100
    vfill "f" using 1:2:3:(0):4

Here we have set the voxel grid ranges (vxrange, etc.) to be the same as the coordinate ranges in the data, and they have the same dimensions, so the data will correspond exactly to the voxel grid. In vfill we use a radius of 0 so that each point in the "source" will only affect the voxel that it lies on top of. In a using clause, bare numbers refer to columns in the data, but parentheses are used to delimit arithmetic expressions, so (0) is just the number zero.

For our dipole example, we want the voxel data to come from the value of the potential, which you may remember looks like charge/r, where r is the distance from the charge (we are not worrying about units here). But that calculation goes infinite when r is zero, so we need to handle that case, which we can do by using gnuplot's ternary operator to define the potential function as:

    pot(r) = r > 0 ? 1/r : 10**6

which avoids the infinity by substituting an arbitrary large number.

We will set the ranges of the bounding box and the voxel grid to go from 0 to 1, and in the vfill command we will use a radius large enough to include the entire box, because the range of the potential field is infinite:

    vfill $charges using 1:2:3:(2):($4*pot(VoxelDistance))

Here we have used gnuplot's convenient VoxelDistance function, which computes r for us. In the parenthesized expressions in using clauses, source columns are referred to with dollar signs; so here we are extracting the charge values from the fourth column in the data. After executing the vfill, we will be able to plot the voxel grid in various ways.

It would make sense to use a color palette for representing the potential field that clearly distinguishes between positive and negative values, and that blends into the white background for values close to zero. We can define a palette with these properties with the command:

    set palette define (0 "red", 0.5 "white", 1 "blue")

This maps the entire range of values stored in the voxels to the range zero to one, then uses those mapped values to determine which color to choose from the palette. So this will map the smallest (most negative) voxel values to pure red, blending into white in the middle of the range, where the potential is zero, and grading to pure blue at the maximum. We have already used this palette to plot the two charges in our first illustration.

Now the command

    splot $v with points above -10**6 pointtype 7 pointsize 0.4 linecolor palette

will plot the voxel grid $v with a style that places a point at each voxel location. The above clause plots all those points with values above the given number; as there is no built-in way to just plot everything, we are obligated to choose a number smaller than any that exists on the grid. The pointtype 7 means closed circles, and the "linecolor palette" tells gnuplot to color the points using the palette that we defined above.

This plot suffers from a kind of moiré pattern, due to viewing periodic arrays of points superimposed on each other. This artifact is not only annoying, but can give a false impression about the nature of the data that the visualization is trying to clarify. Fortunately, the problem can be eliminated by applying jitter. The jitter command appeared first in gnuplot 5.2, and was described here in an earlier article. It applies small random displacements to plotted points; this application of noise is just what is needed to break up the alignments that cause the moiré interference in voxel plots.

If we insert the following command before the splot command in the aforementioned script:

    set jitter spread 0.3

we get the following output:

There are many possible variations of the point-plotting of voxel grids. The opacities, colors, types, and sizes of the plotted points can all be made to depend on the field values in various ways. This can quickly become complicated, but we have one example in the Appendix.

Projection onto surfaces

Rather than trying to visualize the whole volume all at once with a mass of points, a popular alternative technique for 3D visualization is to cut the volume with a plane and plot the projection of the field on it. This is done using the "++" synthetic data source, which will generate coordinate pairs over the bounding box range, with steps based on the samples and isolines settings. It is explained, tersely, on pp. 112–113 of the gnuplot manual [2.1MB PDF]. The command:

    splot "++" using (0.2):1:2

will draw a plane at x = 0.2, parallel to the y-z axis.

We want to color this plane according to the values of the voxel grid elements that it intersects. To do this, we add a fourth field to the using clause that takes the values from the voxels:

    splot "++" using (0.2):1:2:(voxel(0.2, $1, $2)) with pm3d

With the same preliminaries as before, and after making the settings mentioned above, the command produces this plot:

To get an overview of the 3D structure of the field, we need more than one slice. By looping through a series of splot commands, we can lay down a set of slices. If we make the surfaces transparent, we can create a clear visualization of the field. To do this, we'll replace the single splot command above with this:

    set style fill transparent solid 0.4
    splot for [j=0:9] "++" using (j/10.):1:2:(voxel(j/10., $1, $2)) with pm3d

The solid 0.4 sets the opacity of the surfaces; the smaller the number, the more transparent they are. The splot loop replaces the constant x-value we used for the single slice with 10 values from 0 to 9 (thus generating slices at each x, 0.0 to 0.9).

While probably not necessary for this example, the visualization of complex 3D structures can sometimes be aided by animating some aspect of the plot: the angle of view, the range of values plotted, or the position of a cross-section. We can easily adapt the previous command to save a series of plots on disk, using the block iteration syntax introduced in gnuplot 4.6:

    do for [j=0:99] {
        set out "slfr" . gprintf('%02g', j) . '.png'
        splot "++" using (j/99.):1:2:(voxel(j/99., $1, $2)) with pm3d
    }

This uses the same splot command as in the immediately preceding example, but with 100 slices rather than 10. Before each splot command we set the output file name to something that contains the slice number using two digits, so that the names will be sorted correctly for post-processing. We construct the file name string with gnuplot's gprintf() function, which is similar to the familiar sprintf() function from C (and which is also available within gnuplot), but a little simpler. After running this, we will be in possession of 100 files named slfr00.png to slfr99.png.

There are many programs that can stitch these up into a movie. For these purposes I use ImageMagick. If you have this installed, the command is simply:

    convert slfr*.png slfr.gif

And the result is:

Isosurfaces

The splot command for visualizing voxel grids has one more trick up its sleeve: isosurfaces. This is a surface, or a set of surfaces, that show where the voxel grid has certain values. In other words, it is the 3D analogue of contour lines for functions of two variables. And, as with contour lines, they may need to be interpolated between grid values.

The following command draws the isosurfaces showing where the voxel grid has the values 1 or -1. Gnuplot colors the "top" and "bottom" of the surfaces differently by default. The rest of the setup is the same as before, so we need not repeat it—but the depthorder command is necessary for the rendering to be correct; it ensures that pieces of the surface "closer" to the viewer are drawn over those farther away.

    set pm3d depthorder
    splot $v with isosurface level 1, $v with isosurface level -1

That results in the following graph:

We can also apply transparency and iteration to isosurface plotting, to create a set of nested isosurfaces that amount to a 3D contour plot, perhaps the most useful visualization for this particular problem:

    set style fill transparent solid 0.4
    splot for [j=-100:100:1] $v with isosurface level j

Here we see all of the isosurfaces for the whole number values from -100 to 100:

Visualization of 3D data is inherently complicated. Admittedly, some of the commands and concepts that gnuplot exposes for voxel plotting can take a while to understand; that was certainly the case for me. The official manual is quite terse; most of what can be done with voxel grids I discovered through experimentation and brief excursions through the source code. But one reason that these commands seem arcane is their flexibility; once they are mastered they enable almost anything that can be imagined. For example, the technique described here for cutting the space with a plane can be immediately extended to the projection of voxel data onto curved surfaces, such as a cylinder embedded within the potential field. In other words, gnuplot may sometimes demand a bit much from the user, but it offers a lot in return.

Index entries for this article
GuestArticles	Phillips, Lee

---

to post comments