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.

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 2Then 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+1Where '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+1The program fragment on the right implements the useful

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 bisectionOnce 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

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 bisectionsBisections 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:

The equations on the right give the values of a , b , c . Hence the vectorsra =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)

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 IsosurfOutput 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:

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

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.

- Contours running in clockwise direction. They include a
**valley**. - Contours running in counterclockwise direction. They include a
**hill**.

- ==== : NOP
- ===< : coplanar with boundary
- ===> : coplanar with boundary
- ==<< : NOP
- ==<> : through one edge
- ==>> : NOP
- =<<< : NOP
- =<<> : through one vertex
- =<>> : through one vertex
- =>>> : NOP
- <<<< : NOP
- <<<> : triangular surface
- <<>> : triangular surface
- <>>> : quadrilateral surface -> two triangles
- >>>> : NOP

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.