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 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 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.psSo far so good for the Fortran implementation of Paul Bourke's algorithm. Meanwhile, my favorite programming language has become (some good old version of) Delphi Pascal. So here comes my most recent implementation of the algorithm in that language : PaulBourke.pas .
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:
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 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: no care has been taken of the direction of the normals. Hence you should always switch on 'bi-directional' lights.
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.
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.