I am proposing here that we develop our own isosurface generator for use in the cave. The surfaces I have seen generated by Tecplot appear to have too many small patches in regions where curvatures of the surface are relatively small, and larger patches would do. I think we need a surface generator that has the following properties:

• Patches on F(x,y,z)=0 (say) are defined within an unstructured collection of polyhedrons that fill the appropriate parts of 3-dimensional space.
• Values of F, x, y and z should be obtained by interpolation, only in parts of the underlying data sets where F is near zero.
• Subdivision of the polyhedrons to reduce the sizes of patches should be relatively easy to implement.
• Automatic refinement in regions of high curvature should be implemented.
• Algorithms for generating the patches should not have to treat any ambiguous cases.
• It should be possible to use parallel processing in all parts of the package.

What I have to offer at this point is an approach to some of the points outlined above. What we need is a lot of help with the other parts and help with translation of my MATLAB prototypes into more suitable implementations.

Unambiguous decisions follow from two familiar devices:

• Accidental zeros of F that occur on vertices of the polyhedrons are perturbed to arbitrarily small values, positive or negative, randomly or not. This reduces the number of ways a patch near F=0 can intersect a polyhedron.
• The polyhedrons are tetrahedrons (not necessarily regular) and pyramids (the egyptian kind with five faces). These two are the smallest ones for which the disposition of sign changes on edges defines the orientation of a surface patch uniquely, and it is easy to fill space with them. The cases are:
  
  tetrahedron
  
  three zeros -> triangle (4 ways)
  four zeros -> quadrilateral (3 ways)
  
  pyramid
  
  three zeros -> triangle (4 ways)
  four zeros -> quadrilateral (5 ways)
  five zeros -> pentagon (4 ways)
  six zeros -> two triangles (2 ways)

By way of an example, here is a typical set of isosurfaces.

The function is

F(x,y,z) = x2/4 + sin2y - sin2z

and the surfaces are F=-0.5, F=0 and F=0.5.

The surfaces are composed of 3-sided (red), 4-sided (green) and 5-sided (blue) patches. They have been constructed from numerical values of F,x,y and z in simple-cubic arrays. The (MATLAB) prototype of a surface generator of the kind I am suggesting here produced the data for the 18078 patches in 2.6 seconds on a fairly slow (266 Mhz) PC. Similar algorithms that use only tetrahedrons produce triangles and quadrilaterals at a rate slightly greater than 8000 patches per second.

The algorithm that produced the vertices on these patches where F is near zero (say) is one of several different ones that subdivide elementary (empty) cubes that contain a part of the surface into tetrahedrons and/or pyramids (five faces). It proceeds as follows

• The entire array F is first scanned for zeros with »N=find(F==0) (the Matlab prompt » indicates a line of code), and the accidental zeros that are at vertices are perturbed with »F(N)=eps*(2*round(rand(size(N)))-1).
• The elementary cubes are then scanned for changes of the sign of F on the edges. There are none with one or two sign changes, and only those with three or more contribute to the collection of patches.
• In the elementary cubes that contain a part of F=0 the averages of x, y, z and F are used to construct data at a central node for a set of six pyramids that fill the cube. (Variants of this use tetrahedrons or both tetrahedrons and pyramids.)
• The pyramids that may contain a part of F=0 are scanned for accidental zeros (which are perturbed) at the apex, and the result is a collection of pyramids that can have 0, 3, 4, 5 or 6 sign changes of F on the edges.
• The utility in pyras.m encodes the 6+5+4 ways there can be patches with 3, 4 or 5 sides and does linear interpolations to find the vertices of the patches.

The point of isolating the computations in pyras() and a similar utility tetras() was to have them all done in an unstructured context. The cubic array of data points in the example was incidental. Advantages of this are:

• Data from any collection of pyramids and/or tetrahedrons can be used by pyras() and/or tetras() to find patches near F=0.
• There are fairly convenient subdivisions (refinements) of the tetrahedron and the pyramid, some with tetrahedrons and pyramids and others with tetrahedrons alone.
• The encodings of cases and construction of patches is strictly deterministic when accidental zeros have been perturbed.
• The code is relatively efficient, and inherently parallel.
• The result is altogether independent of the topology of the surfaces. (Of course that could be a disadvantage too.)

Here again is a zoomed view of F=[-0.5 0.5] in which the isosurfaces have been partially, doubly covered and sierpinskied.

Here is a slightly different view in which pentagons have been left out and triangles and quadrilaterals have been shrunk.

And here it is as Roy Lichtenstein might have done it.

Or maybe Georges Seurat.

Here are some of the simpler refinements:

The truncated tetrahedron is an octahedron. The three planes cut it in two pyramids, if used one at a time, or in eight tetrahedrons. Another possibility is to choose one of the three diagonals to be the common edge of four tetrahedrons. The eight tetrahedrons then fill equal volumes in the computational grid.

Easy ways to slice a pyramid use a diagonal of the base to define two tetrahedrons or the centroid of the base to define four. There are many other possibilities, including this, with six pyramids and four tetrahedrons.

When the edges of the pyramids have equal lengths, as shown here, and the refinements are tetrahedron -> four tetrahedrons & two pyramids and pyramid -> six pyramids & four tetrahedrons, the tetrahedrons are regular and the pyramids are all geometrically similar.