sci.math.numanalysis
SUNA55, Contours & Isosurfaces
==============================
A nice article, describing a robust algorithm for calculating the contours of a
function z = F(x,y) , is written by Paul D. Bourke: "A Contouring Subroutine",
BYTE journal, June 1987. Generalization of this concept of contouring to three
dimensions is known to be the construction of "isosurfaces". A leading article
describing a procedure for doing this is: William E. Lorenzen, Harvey E. Cline,
"Marching Cubes: a high resolution 3D surface construction algorithm", Computer
Graphics, volume 21, Number 4, July 1987. Based upon the theory in this article
is the "Marching Cubes surface tiler" I once received from Steve Lamont. Quoted
with his permission (: I'm pretty sure):
> North Carolina Supercomputing Center
>
> March 1990
> The output from mcube looks like
>
> xmin xmax ymin ymax zmin zmax
> x1 y1 z1
> x2 y2 z2
> x3 y3 z3
> x4 y4 z4
> 
> ...
> x1 y1 z1
> x2 y2 z2
> x3 y3 z3
> 
>
> where the [xyz]min and [xyz]max are copied from the input (and ignored by
> the program unless the bounding box option is chosen  see below for usage)
> The xn, yn, zn values are vertex coordinates. In the case of quadrilaterals
> or higher order polygons, the vertices are *not* necessarily coplanar and must
> be post processed by a triangulator filter if planar polygons are required.
> Each polygon is terminated by a single '' (minus sign) followed by a linefeed
We will call this "mcube" output format "spl". Files with this format will have
".spl" suffix. They are generated by our own surface tiler (: see below).
I guess that also thee module "isosurface" in AVS (: Application Visualization
System) is 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 have problems with the latter thing. At least with ConvexAVS 3.9 it
gives results which cannot possibly be correct. That's one reason why I decided
to devise my own theory & code. However, the approach in this SUNA 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.
The basic element for contouring a la Bourke is a triangle. The function value
at an isoline can be set to zero without loss of gererality. With respect to
this isovalue 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 in numbers.
A correspondence between the inequalities 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
be written in Fortran 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=more1  m=0
do 15 k=1,3  do 15 k=1,3
num(4k)=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 function 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", but I (well, almost) forgot
how it works;). Anyway, the output of "permut" is stored in just another array,
called np(3,6) . Once we have found the ternary decomposition of an index, all
permutations of it are given by num(np(j,i)) , and indices "more" belonging to
these permutations are 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
whichareessentiallydifferent are set to zero.
Loop through all indices L from 1 to 27. If must(L) or move(L) is already
nonzero, 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 we determine later on
the function values fv(j)=f(np(j,i)) , they must be transformed _back_ to the
case L we started with. The inverse of an inverse is the thing itself, right?
The inverse "no" of a permutation "nr" is found by no(nr(k))=k , for k=1,..,3 ,
almost by definition. So far so good.
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(3num(k),3) .
The above procedure has been implemented in subroutine "setup3". It results in
six cases which are essentially different, and for which a separate calculation
of an isoline has to be devised. The "setup3" routine is called only _once_ &
it 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, eh? ;).
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" differs little from that.
Quadrilaterals can be built up from triangles in several ways. I have some fun
in doing it here in a rather deviant way, namely four _overlapping_ triangles.
This results in "not so nice" contour plots, but I like them: because you can
see very clearly the inaccuracies involved. No faking.
Oh yeah, the accompanying computer program is written in Fortran, it's available
at the well known "ftp" site, and it's called "suna55a.for", together with some
sample data "suna55a.dat" (which, by the way, is the streamfunction, calculated
for ideal flow past a circular cylinder).
Let's go to the 3D case now: isosurfaces. 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 even difficult to generalize this administration to multiple dimensions .. !
"setup4" results in NINE cases which are essentially different, and for which a
separate calculation procedure of an isosurface has to be devised. The essential
different cases are:
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 sides 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 left give the values of a , b , c . Hence the vectors r
pointing to the points of intersection can be calculated.
The kind of calculations that are contained in our subroutine "tetrahedron".
Subroutine "tetrahedron" is called in turn by "brick", which makes isosurfaces
in a brick (indeed), by subdividing it into five (nonoverlapping) tetrahedra.
I know that more advanced subdivisions would be possible, for example 24 pieces
by employing midsurface and midbulk values. "Brick" also works for degenerate
cases, such as prisms. A routine called "volume" takes care of volumes equal 0 .
Subroutine "brick" is embedded in the routine that organizes the bulk of our
calculations: "bereken" (: Dutch for "calculate"). Which at last is called by
the main program. The file "suna55b.for" contains everything.
Files in AVS UCD format, detectable by their suffix "inp", are the input files
that are understood by "suna55b.for". Output is in "Steve P. Lamont" format,
which was described above. These "spl" files can be converted into AVS "Geom"
files by a short program, which is added as "suna55c.for" to the collection.
The example file "suna56a.inp" represents a discretization of the potential of
ideal flow past a sphere, publication of which is forthcoming as "suna56.net".

* Han de Bruijn; Applications&Graphics  "A little bit of Physics * No
* TUD Computing Centre; P.O. Box 354  would be NO idleness in * Oil
* 2600 AJ Delft; The Netherlands.  Mathematics" (HdB). * for
* Email: Han.deBruijn@RC.TUDelft.NL  Fax: +31 15 78 37 87 * Blood