Contours & Isosurfaces

A robust algorithm for calculating the contours of a function z = F(x,y) has been published by Paul Bourke. Generalization of the concept of contouring to three dimensions is known as the construction of 'iso-surfaces'. The leading algorithm for doing this is 'Marching cubes'. Based upon this is the surface tiler I once received, for free, from Steve Lamont. (Thanks Steve!)

My guess is that the module 'isosurface' in (for example) AVS is also based upon the Marching Cubes algorithm. The UCD counterpart of this module, which is meant for Finite Elements instead of Finite Volumes, is called 'ucd iso'. I had problems with the latter thing. At least with ConvexAVS 3.9 it gave results which cannot possibly be correct. That's one reason why I decided to devise my own theory & code. However, the approach in this article will be different from Marching Cubes. It is based instead upon a generalization of Bourke's 2D contouring method. So let's take a close look at contours first.

Iso-lines

The basic element for contouring a la Bourke is a triangle. The function value at an iso-line can be set to zero without loss of generality. With respect to this iso-value 0 , the function values at the three vertices, numbered 1,2,3 , can be equal, less than, or greater than zero. For example let f1 = 0 , f2 < 0 and f3 > 0 . We abbreviate this by '=<>' . Likewise '<<<' for f1, f2, f3 < 0 . Doing this systematically, we get the following, exhausting, list of 27 cases:
      1    2    3    4    5    6    7    8    9
     ===  ==<  ==>  =<=  =<<  =<>  =>=  =><  =>>

     10   11   12   13   14   15   16   17   18
     <==  <=<  <=>  <<=  <<<  <<>  <>=  <><  <>>

     19   20   21   22   23   24   25   26   27
     >==  >=<  >=>  ><=  ><<  ><>  >>=  >><  >>>
The above already represents an encoding of logical decisions into numbers. A correspondence between (in)equalities and digits is established as follows:
     replace  =  by  0
     replace  <  by  1
     replace  >  by  2
Then for example =<> corresponds with 012 , which is a number coded base 3 : 012 = 1*3 + 2 = 5 decimal. Adding 1 to this result gives 6 , which indeed is the index of =<> in the above table. The corresponding program fragment could have been written, in Fortran for example, as follows:
      more=0
      do 10 k=1,3
      more=more*3
*     if(f(k).eq.0.) more=more+0
      if(f(k).lt.0.) more=more+1
      if(f(k).gt.0.) more=more+2
10    continue
      more=more+1
Where 'more' now serves as an index for further actions to be taken.

Let's consider the (base 3) numbers: 012 , 021 , 201 , 102 , 120 , 210 . They correspond with the following cases: =<> , =>< , >=< , <=> , <>= , ><= . It is clear that these cases are essentially 'the same', to be precise: apart from a permutation. Thus we may build all the permutations of the digits in a number 012 , and identify these numbers in some way. Finding the base 3 digits in all indexes 'more', as defined above, is done by the program fragment on the left:

       m=more-1                    |       m=0
       do 15 k=1,3                 |       do 15 k=1,3
       num(4-k)=mod(m,3)           | 15    m=m*3+num(k)
 15    m=m/3                       |       more=m+1
The program fragment on the right implements the useful inverse of it.

The task is now to find all permutations of an array num(3) . I have a little routine for doing just that. It's called 'permut' (one of my beloved routines, but I have almost forgotten how it works ... ;-). Anyway, the output of 'permut' is stored in just another array, called np(3,6) . Once we have found the (ternary) digits of a given index, all permutations of it are given by num(np(j,i)) , and indices 'more' belonging to these permutations will be constructed by the above inverse function.

Now the overall procedure becomes clear. At first, define two arrays: must(27) for all cases which are essentially different, and move(27) indicating which permutation has to be used. Initially, both arrays and a counter for the cases which-are-essentially-different are set to zero.

Loop through all indices L from 1 to 27. If must(L) or move(L) is already non-zero, do nothing. If not, increment the 'essential case' counter, and loop through all permutations of the (ternary) digits of L . Calculate accompanying indices 'more' and store the value of the case counter in must(more) . Also, store the index of the permutation in move(more) . Huh? Store the index of the inverse of the permutation in move(more) ! Because, if the function values fv(j)=f(np(j,i)) are to be determined later on, they must be transformed back to the case L we started with. The inverse of an inverse is the thing itself, right? Now it should be well known that the inverse 'no' of a permutation 'nr' is found by no(nr(k))=k , for k=1,..,3 , almost by definition.

There is another class of cases to be identified, however. According to both the articles of Lorenzen/Cline and Bourke, it makes no difference if we exchange the inequalities < and > , or our digits 1 and 2 . But digits 0 should be preserved, as in: num(k)=mod(3-num(k),3) .

The above procedure has been implemented in a subroutine ' setup3'. It results in six cases which are essentially different, and for which a separate calculation of the iso-line has to be devised. This 'setup3' routine is called only once and could entirely be replaced by data statements for the tables that are set up by it (but then a nice piece of theory would probably be lost, huh? ;-). OK. The essential cases are listed here:

label   logic   base 3   decimal   Bourke   isoline
-----   -----   ------   -------   ------   -------
   1.    ===     000        1        VI      NOP
   2.    ==<     001        2         I      draw from vertex to vertex
   3.    =<<     011        5         V      NOP
   4.    =<>     012        6        II      draw from vertex to bisection
   5.    <<<     111       14        IV      NOP
   6.    <<>     112       15       III      draw from bisection to bisection
Once we have setup the administration properly, the rest of the coding becomes simple & straightforward. It's worked out further in Bourke's article, and the implementation in our subroutine 'triangle' is not much different from that, except for the fact that, thanks to the above, it's kindof 'IF-less' programming, with 'cases' (Fortran: computed goto's) instead.

PostScript will be selected as the preferred graphical output format. A couple of routines are associated with this choice, namely for initialization and for drawing. Limiting the number of characters in the 'contour.ps' file is accomplished by a base 10 logarithm function. The plotfile is connected to Fortran unit 7 .

If nothing else is done, then a lot of senseless computation will be involved when invoking subroutine 'triangle'. The problem is that most of the contour levels are simply not present at a single triangle. Therefore a routine called bounds was developed, which determines the available range of the contour levels by binary search.

A sample input for the main program can be generated by using a random number generator, followed by a call to 'voronoi'. After that, the contours are generated and (pre)viewed as follows:

 fc Random.FOR
 a.out
 # Outputs: 't' , 'coor2d.dat' , 'coor3d.dat'
 voronoi -t < t > topol2d.dat
 fc Contour.FOR bounds.FOR triangle.FOR setup3.FOR permut.FOR \
 plots.FOR plot.FOR log10.FOR log2.FOR -o Contour
 Contour
 gs contour.ps

Iso-surfaces

Let's go to the 3-D case now: iso-surfaces. Bourke's triangles are replaced by tetrahedrons, with 4 vertices instead of 3 , and 24 permutations instead of 6 .

The enumeration of all cases for isosurfaces in a tetrahedron can easily be done by copying three times the table for the triangle, adding in front of each copy respectively = , < , > , and continuing the numbering. This results in a total of 81 not essentially different cases. The subroutine 'setup3', used for the contouring administration, is copied to 'setup4'. The only changes that have to be made are in the dimensioning of array's and ranges of variables in the loops. Make a 'diff' of 'setup3' and 'setup4' and you will see what I mean. It wouldn't be so difficult to generalize this administration even to ... multiple dimensions greater than 3 !

Invoking the routine 'setup4' results in nine cases which are essentially different, and for which a separate calculation procedure of an isosurface has to be devised. The truly different cases are enumerated below:

label   logic   base 3   decimal   isosurface
-----   -----   ------   -------   ----------
   1.   ====     0000        1     NOP
   2.   ===<     0001        2     triangle coplanar with tetrahedron boundary
   3.   ==<<     0011        5     NOP
   4.   ==<>     0012        6     triangle through 2 vertices and 1 bisection
   5.   =<<<     0111       14     NOP
   6.   =<<>     0112       15     triangle through 1 vertex and 2 bisections
   7.   <<<<     1111       41     NOP
   8.   <<<>     1112       42     triangle through 3 bisections
   9.   <<>>     1122       45     quadrilateral through 4 bisections
Bisections are determined as follows. Take the case with label 8 as an example. The isosurface triangle has it's corners a , b , c at the edges that join the vertices (1,4) , (2,4) , (3,4) of the tetrahedron. The points of intersection are defined by:
ra = r4 + a.(r1 - r4)   with   fa = 0 = f4 + a.(f1 - f4)
rb = r4 + b.(r2 - r4)   with   fb = 0 = f4 + b.(f2 - f4)
rc = r4 + c.(r3 - r4)   with   fc = 0 = f4 + c.(f3 - f4)
The equations on the right give the values of a , b , c . Hence the vectors r pointing to the points of intersection can be calculated. This is implemented in a subroutine called tetrahedron.

What's needed for all routines to work is a three-dimensional mesh, consisting of tetrahedra. As in the 2-D case, such a mesh can be generated by using a 3-D Delaunay triangulator. Use has been made of the previously generated files 'coor2d.dat' and 'coor3d.dat':

 tetras < coor2d.dat > topol3d.dat
 fc Isosurf.FOR tetraedr.FOR setup4.FOR permut.FOR -o Isosurf
 Isosurf
Output of the 'Isosurf' program is in 'Steve P. Lamont' format, which is described elswhere. These 'spl' files can be converted into AVS 'Geom' files by a short program. Note: no care has been taken of the direction of the normals. Hence you should always switch on 'bi-directional' lights.
However, the question whether an arbitrary isosurface is rendered correctly is difficult to judge by visual inspection. Another thing is that a surface tiler for bricks, instead of tetrahedra, is needed for many applications. Therefore an isosurface tiler for a regular surface, namely a sphere, using bricks in a regular grid, is supplied in addition. The example can be executed as follows:
 fc Bol.FOR brick.FOR tetraedr.FOR setup4.FOR permut.FOR -o Bol
 Bol
 metavs bolsurf > bolsurf.geom

Programming

O yeah, the programs contain a number of Fortran 'parameter' statements, by which the size of arrays and the maximum length of files is determined. Their values in the delivery, here, are quite arbitrary: don't forget to adjust them to the dimensions of your actual problem. It can give rise to nasty bugs if you don't!
  Random.FOR:   parameter (LIMIT=20)
 Contour.FOR:   parameter (NODES=1500,LINES=4500)
 Isosurf.FOR:   parameter (NODES=25,LINES=75)
     Bol.FOR:   parameter (NODES=25)
  metavs.FOR:   parameter (LIMIT=100000)
At last, all Fortran programs in this WebPage, and others, have been zipped into one (16 KB) file for convenience.

Pascal / Delphi Updates

Throwing away Directional Information does not much harm with contouring in the 2-D case, at first sight. But if that information is retained, then a quite useful distinction can be made between two kinds of closed contours:
  1. Contours running in clockwise direction. They include a valley.
  2. Contours running in counterclockwise direction. They include a hill.
The distinction is even more crucial in 3-D, with Isosurfaces. One shouldn't feel the need for bi-directional light tricks at all, once the surfaces are rendered correctly, namely: with their normals pointing in the direction of the gradient of the function f(x,y,z). Instead of nine cases, fifteen cases (Equivalence Classes) are distinguished then:
  1. ==== : NOP
  2. ===< : coplanar with boundary
  3. ===> : coplanar with boundary
  4. ==<< : NOP
  5. ==<> : through one edge
  6. ==>> : NOP
  7. =<<< : NOP
  8. =<<> : through one vertex
  9. =<>> : through one vertex
  10. =>>> : NOP
  11. <<<< : NOP
  12. <<<> : triangular surface
  13. <<>> : triangular surface
  14. <>>> : quadrilateral surface -> two triangles
  15. >>>> : NOP
The surfaces are thus oriented in such a manner that their normal is pointing away from the function values '<' and pointing towards the function values '>'.
As has been announced earlier, I am in the process of translating all of my Fortran programs into Delphi (PAscal) 3,5,6 . The Isosurface algorithm is implemented in Turbo Pascal / Delphi 5 nowadays. It can be executed by linking the Unit to a simple TP program:
        program Testen;
        Uses IsoSurf;
        begin
          Kubus;
        end.
A file called 'isocube.spl' is then produced, which in turn can be fed into a renderer which is based upon the Silicon Graphics Open GL runtime library (under Windows NT & Windows 98). The rendering device has been programmed by Stefan. Another useful link is this one.