Searching for vertices of Waterman Polyhedron
As we said already there are known algorithms to produce points in 3D that can be used to create a Waterman polyhedron. However, most of them are to slow to be implemented in MuPAD. For example, the algorithm used in [1] to create Waterman spheres, after transforming it to MuPAD, turned into code with runtime complexity equal to
Such algorithm will work considerably well in binary, i.e. complied, programs. However, programs in MuPAD, like in any other Computer Algebra System, are interpreted and therefore they are much slower. Therefore, we needed a fast and very efficient algorithm. For this purpose it was created an algorithm that cuts a sphere with a given radius R into slices and chooses from them only those points that have coordinates x+y+z divisible by 2. With some improvements done by Oliver Kluge from SciFace Software GmbH. the algorithm produces results in a considerably faster time. For instance in MuPAD in order to produce w1,100, using the algorithm from [1], we need 1,096,586 ms that is about 18 min. While the optimized algorithm produces the same polyhedron in 2,834 ms, that is about 3 sec.
Below, I show the optimized algorithm implemented as a procedure WatermanData.
WatermanData := proc(center, radius:Type::Positive)
local counter, x, y, z, a, b, c, xra, xrb, yra, yrb,
zra, zrb,R, Ry, s, radius2, WD;
begin
a := center[1]:
b := center[2]:
c := center[3]:
WD := table():
counter := 0:
s := radius;
radius2 := radius^2:
xra := ceil(a-s);
xrb := floor(a+s);
for x from xra to xrb do
R := radius2 - (x-a)*(x-a);
if R < 0 then next end_if;
s := sqrt(R);
yra := ceil(b-s);
yrb := floor(b+s);
for y from yra to yrb do
Ry := R - (y-b)*(y-b);
if Ry < 0 then next end_if; //case Ry < 0
if Ry = 0 and is(c, Type::Integer) then //case Ry=0
if (x+y+c) mod 2 <> 0 then
next
else
zra := c;
zrb := c;
end_if
else // case Ry > 0
s := sqrt(Ry);
zra := ceil(c-s);
zrb := floor(c+s);
if (x+y) mod 2 = 0 then // (x+y)mod2=0
if zra mod 2 <> 0 then
if zra <= c then
zra := zra + 1
else
zra := zra - 1
end_if
end_if
else // (x+y) mod 2 <> 0
if zra mod 2 = 0 then
if zra <= c then
zra := zra + 1
else
zra := zra - 1
end_if
end_if
end_if
end_if;
for z from zra to zrb step 2 do
counter := counter + 1:
WD[counter] := [x, y, z]
end_for:
end_for:
end_for:
return( WD )
end_proc:
Observe, the final output from this procedure is a table of point coordinates in 3D. MuPAD tables are very flexible data structures that do not need declarations of size while initiating them.
Using external program to produce convex hull
Producing a table of points was the very first step towards producing a Waterman polyhedron. In the next step we need to filter out those points that do not lie on the surface of the polyhedron and use the remaining points to create convex hull, i.e. the smallest polyhedron containing our points. Instead of implementing a very complex algorithm to produce convex hull, that certainly would be much slower in MuPAD than in binary code, I decided to use a ready program qconvex.exe written in C++. As I did mention before, the QConvex program is part of the Qhull library developed at the National Science and Technology Research Center for Computation and Visualization of Geometric Structures, The Geometry Center, University of Minnesota. The complete Qhull library can be downloaded from http://www.qhull.org.
The simplest way to use QConvex on data produced by the WatermanData procedure is to convert these data to a format acceptable by QConvex, save the data into a text file, apply QConvex on this file and produced results read into MuPAD again.
Converting data produced by WatermanData procedure is quite straightforward task and saving these data into a file is also simple job. In order to call QConvex from inside of MuPAD the shell module was used. Below I show the procedure useQconvex. Observe, the procedure uses a global variable QCONVEX_PATH pointing into a place where the qconvex.exe file was saved in our system, then builds a complete access path to this file, and finally calls qconvex.exe with a list of parameters.
useQconvex := proc() local exe_file; begin shell::changeDir(tempDir): exe_file := QCONVEX_PATH . "qconvex": exe_file := "\"" . exe_file . "\"": shell::system(exe_file . " o <data.txt >data2.txt"): end_proc:
The resulting file data2.txt contains description of the convex hull in Qconvex format. Below I show an extract from such file with my comments.
3 // dimension of the space 201 14 36 // number of points, nr of faces, number of edges -4 -2 0 // here start declarations points -4 -1 -1 -4 -1 1 -4 0 -2 ..... 4 2 0 // end of declarations of points // declarations of faces starts here, first number is the number of points, // used to create a face, remaining numbers are consecutive indexes of points. // For example point [-4,-2,0] has index 0, [-4,-1,-1] index 1, etc. 4 3 8 5 0 4 165 93 39 111 6 5 8 49 118 111 39 6 3 35 107 116 49 8 4 116 175 118 49 6 197 192 151 84 93 165 6 39 93 84 25 0 5 6 35 3 0 25 82 89 4 82 25 84 151 6 192 195 161 89 82 151 4 89 161 107 35 6 118 175 200 197 165 111 6 107 161 195 200 175 116 4 200 195 192 197
The next step would be to use these data to produce faces of a polyhedron.
The SurfaceSet class and its role in Waterman polyhedra creation
It is obvious that faces of a polyhedron are polygons in 3D. We can easily create them using the SurfaceSet plot class. Although the SurfaceSet class has many different options, we needed only the one that takes coordinates of points and builds a sequence of triangles combining these triangles into a polygon. If these points lay on a common plane then we will obtain a flat face of a polyhedron. Here is a very simple example showing how SurfaceSet class works. I deliberately use here points that do not lie on a common plane, and I also show edges of triangles. This will give us more information on how SurfaceSet is created.
data := [
1.5, 2.6, -1.0,
-1.5, 2.6, 1.0,
-3.0, 0.0, -1.0,
-1.5, -2.6, 1.0,
1.5, -2.6, -1.0,
3.0, 0.0, 1.0
]:
plot(
plot::SurfaceSet(data,
MeshListType=TriangleFan,
MeshVisible=TRUE,
Axes=Boxed
)
)
Observe, by default SurfaceSet will not produce edges of triangles, as well as edges of the polygon. Therefore, if we wish to see edges of a polyhedron we have to create them separately as an independent plot object.
In the final creation we produce two groups of objects: faces and edges of the polyhedron. These two groups form a plot object that can be displayed on a computer screen using the plot command.
One of the special features of this implementation are random flat colors. The major reason for this approach was to avoid the general MuPAD coloring scheme with red color in top parts of the graph, and blue color in the bottom. MuPAD default coloring scheme works well with smooth surfaces, but at the same time makes very boring images of polyhedra. Here is the code showing how random colors of faces were created.
face := plot::SurfaceSet(facePoints, MeshListType=TriangleFan, FillColor=[frandom(),frandom(),frandom()], FillColorType=Flat ):
Of course, this place is a great opportunity to experiment different coloring schemes. We may return to this problem in another article.
Developing Waterman polyhedra was a very interesting computational problem. The methods used here can be applied to many other similar problems. For example, one may consider completely different set of points in 3D as a starting point for developing polyhedra. Instead of taking a big sphere with radius R, one could take a large octahedron or dodecahedron. This would perhaps lead to a new family of polyhedra. Finally, using external executable programs can enrich calculations in MuPAD and allow access into areas where MuPAD programs doesn't fit or simply were not yet developed.
Literature
- Bourke P., http://astronomy.swin.edu.au/~pbourke/geometry/waterman/
- Newbold M., http://dogfeathers.com/java/ccppoly.html
- Urner K., http://www.4dsolutions.net/ocn/
- Waterman S., http://www.watermanpolyhedron.com/
- Waterman S., Waterman Polyhedra, Symmetry: Culture and Science, Vol. 13, No.1, 2002
- Webb R., http://web.aanet.com.au/robertw/Stella.html
