sci.math.num-analysis 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 "iso-surfaces". 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 iso-line 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=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 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- 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 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(3-num(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 iso-line 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 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 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 (non-overlapping) tetrahedra. I know that more advanced subdivisions would be possible, for example 24 pieces by employing mid-surface and mid-bulk 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 * E-mail: Han.deBruijn@RC.TUDelft.NL --| Fax: +31 15 78 37 87 ----* Blood