O. Michael Melko

Rendering Solid Models with the Aid of 3D Printers

3D printers are a potentially useful tool for geometric visualization in mathematical research and education. In this article, we describe the mathematics of minimal surfaces in some detail, and then we present a Mathematica package for generating solid model data of such surfaces. In particular, we show how the package was used to generate a 3D model of Costa’s surface.


In [1], Palais describes a program for the visualization of mathematics through the use of computer graphics. He notes that visualization has been instrumental in some important mathematical discoveries and is also useful for educational purposes. With this in mind, he proposes the creation of an online interactive gallery of mathematical visualization, which he calls a “mathematical exploratorium.” As helpful as computer graphics are for visualization, it can be argued that there is an additional benefit to be had in viewing and handling an actual physical object. Indeed, there has been a long-standing tradition at German universities of producing plaster models of interesting geometric objects. Nowadays, physical models can be easily created by means of stereolithography, or “3D printing technology.”

3D printers produce solid objects from appropriate input data. They were originally created for rapid prototyping of new product designs but are increasingly being used for other purposes, such as highly customized manufacturing and scientific visualization. Their applications will continue to grow as the underlying technology improves and decreases in cost. One type of 3D printer uses a powder-binder technology to create objects via a layering technique: a thin layer of powder is spread across a planar surface, and then a print head applies a binder within the cross-sectional area of the object being created. This process is repeated, adding layer upon layer, until the object is complete.

As input, 3D printers require data specifying the vertices, polygons, and normals of the object to be rendered. If the object is to be colored, the colors of the polygons must also be provided. Various file formats may be used to store this data, including the Polygon File Format (or PLY), which is also known as the Stanford Triangle Format. This is the file format used to render Costa’s surface and is described later.

We generate point (or vertex) data for Costa’s surface by means of its Weierstrass representation. Loosely speaking, the Weierstrass representation provides a recipe for creating a parametrization of a minimal surface in from two meromorphic functions defined on a Riemann surface. In particular, when these functions are and , where is the Weierstrass function and is a certain constant, the underlying Riemann surface is a torus of the form , where is a discrete lattice in the complex plane . The resulting minimal surface has often been rendered in 3D graphics images. To produce the coordinate data for this surface, we must integrate certain rational functions of and . This results in coordinate functions that are expressed in terms of the Weierstrass and functions. Since these functions are built into Mathematica, it is easy to generate the required data.

To generate a solid model, we must produce vertex and face data for a polyhedron that bounds a volume in . To achieve this, we first choose a proper subregion of a fundamental domain of the functions and that contains no poles. This ensures that the corresponding piece of Costa’s surface is of finite extent. We then generate two surfaces by means of normal displacement and “glue” the resulting boundaries together. This data is then exported to a PLY file, which is used to print the model.

In what follows, we first provide some background from minimal surface theory. This includes a review of classical surface theory, a description of the Weierstrass representation, a summary of pertinent facts about elliptic functions, and a description of the parametrization of Costa’s surface that we use to generate model data. This is followed by a description of the Minimal Surfaces package developed to render the model data. Then we illustrate how to use the functionality provided by this package to create both graphics objects and PLY data. Finally, we discuss some ways in which the work in this paper might be extended, including ideas for mathematical experimentation and enhancements to the Minimal Surfaces package.

A Thumbnail Sketch of Minimal Surface Theory

Intuitively, a smooth surface in Euclidean space is locally area minimizing if any small deformation of results in a surface of larger area. (The precise mathematical definition of minimal surface requires the introduction of some technical preliminaries and is given later.) Soap films spanning a curve in , for example, satisfy this property. In general, may have self-intersections, in which case the surface is said to be immersed, otherwise, we say that it is embedded. The purpose of this section is to describe the Weierstrass representation for minimal surfaces, which provides a way of explicitly constructing a large class of such surfaces via complex function theory. As an example, we shall discuss Costa’s minimal surface in some detail. We begin by summarizing essential facts from classical surface theory; further details may be found in [2].

Essential Facts from Classical Surface Theory

We parametrize a smooth surface in by means of a map , where is an open connected subset (or region) in , and we assume that . We refer to as the trace of . Furthermore, we use to denote a system of coordinates on , and we assume that is smooth, i.e., differentiable to arbitrary order with respect to and . We also assume that is regular, i.e., that the tangent vectors and are linearly independent for every point .

Let denote the standard inner product on , and fix an orientation. Then the metric, or first fundamental form, on with respect to is given by the symmetric tensor

where the coefficients , , and are functions on given by

Let denote the sphere of unit radius in centered at the origin, and let be the unit normal vector at the point in that is consistent with the chosen orientation of . We use to denote the unique element of that is parallel to . The map is called the Gauss map—it provides a measure of how the surface bends in its ambient space. The Gauss map is used to define the second fundamental form,

where the coefficients , , and are again functions on given by

We identify the first and second fundamental forms with the symmetric matrices that define them:

The forms and encode intrinsic and extrinsic geometric properties of the surface . Intrinsic properties are those that are derived from the presence of a distance measure on and do not change under isometric (or distance-preserving) deformations in the ambient space . Extrinsic properties are those that depend on how the surface is immersed into . (As a motivating example, consider wrapping a geographical map into a tube: distances between points within the map do not change, but the way they lie in space does.) A fundamental fact from classical surface theory is that, up to rigid motion, the forms and completely determine the geometry of and how it lies in .

The mean curvature and the Gauss curvature are defined by


Both and are independent of the choice of coordinates on , and it turns out that is an extrinsic property of , while is an intrinsic property. The latter fact is the well-known Theorem Egregrium of Gauss.

An important example of an intrinsic property of is its area, which we denote by . This area is expressed in terms of the first fundamental form as follows:


where is the infinitesimal element of area on with respect to the parametrization . Note that the integral in equation (2) is independent of the choice of parametrization, so that it only depends on the trace of .

Suppose now that is a smooth function. Then we can use to define a normal variation of as follows:

For small values of , the image is a smooth surface near . If denotes the area of , then a straightforward calculation shows that

where and denote the mean curvature and element of area of . If a surface is locally area minimizing, we expect the derivative to vanish for all choices of . This can only happen if vanishes identically. Hence, we have the following:

Theorem. If the surface is locally area minimizing, then the mean curvature of vanishes identically ().

The standard definition of minimal surface is motivated by this fact:

Definition: The surface is minimal if its mean curvature vanishes identically.

Thus is a necessary, but not sufficient, condition for to be locally area minimizing. Determining whether a minimal surface is actually locally area minimizing would entail calculating the second variation of the area functional on . The implication would follow if the second variation were always positive.

The total curvature of is defined to be

We say that has finite total curvature if . Furthermore, we say that is complete if all its geodesics can be extended indefinitely. A surface is said to be of finite topological type if it can be smoothly deformed into a compact surface of finite genus, possibly with several holes. It was long conjectured that the only complete, embedded minimal surfaces in of finite topological type are the plane, the catenoid, and the helicoid. Costa’s minimal surface was the first counterexample to this conjecture to have been found (see [3] for details).

The Weierstrass Representation for Minimal Surfaces

It turns out that the geometry of minimal surfaces is intimately related to complex function theory. This connection leads to a simple recipe for constructing minimal surfaces, which we describe here. We only state the necessary results; further details may be found in [2].

We identify with the complex plane by means of the usual correspondence . Suppose that is a simply connected region in , that is, a region in in which all closed curves can be contracted to a point. A complex-valued function on is said to be holomorphic if its complex derivative exists for all .

Theorem. Suppose that are two holomorphic functions on a simply connected region , and define to be the holomorphic curve with components

Then we have the following:

(i) Componentwise integration


yields a holomorphic curve .

(ii) For each , the trace of the map


is a minimal surface. The collection of all such maps is called the associate family of .

(iii) For any , is the stereographic projection into of the Gauss map of .

Since and are assumed to be holomorphic on , and is assumed to be simply connected, it follows from basic complex function theory that the integration in equation (4) is path independent.

We refer to a triple satisfying the above conditions as the Weierstrass data for the corresponding associate family , and we refer to as the Weierstrass representation of the minimal surface . Note that is an isothermal parametrization for each , that is, the coefficients of the first fundamental form satisfy and In fact, it can be shown that


where . Also, we have


Here, denotes the complex norm of and .

Meromorphic Functions on Complex Tori

Our goal in this subsection is to introduce the Weierstrass data used to obtain a parametrization of Costa’s surface. Before doing so, we provide a little background in elliptic function theory. Details may be found in [4].

Suppose that and are two complex numbers such that . Then the lattice of points is a subgroup of the group of translations on and is isomorphic to the additive group of Gaussian integers . Thus, acts on by the rule , and the quotient space is topologically a torus, which inherits a complex structure from . Let denote the corresponding projection map. We refer to and as basic periods of and sometimes write for .

Any function on lifts to a function on . Such a function satisfies for all and is said to be periodic. It is not possible for a complex function to be both periodic and holomorphic in all of , but there is a rich theory of functions that are periodic and meromorphic. A function is meromorphic on if it is holomorphic on , where is a discrete set without accumulation points in , and if, for any , has a pole at , that is, has a power series expansion in a neighborhood of of the form

The positive integer is called the order of the pole at . A pole is said to be simple if it is of order one. Functions that are both meromorphic and periodic with respect to some lattice are referred to as elliptic functions.

Our objective in what follows is to describe a particular solution of the period problem. In the context of complex tori, it may stated as follows:

The Period Problem: Find meromorphic functions and , periodic with respect to some lattice in , such that is also periodic, where is given by equation (4). Here, the Weierstrass data is , where is a fundamental domain of in , and is the set of poles of and in .

Note that, by making appropriate cuts in , we may consider it to be simply connected. Any solution to this problem will topologically be a torus (possibly with several holes) immersed in .

Costa’s surface arises from what is arguably the most basic of elliptic functions: the Weierstrass function. It is defined by the series expansion


where denotes the set of nonzero elements in . This function has poles of order two at each of the lattice points in and is holomorphic everywhere else in . It therefore projects to a meromorphic function with exactly one pole of order two in . It is known that any meromorphic function on may be expressed as a rational function of and its complex derivative . In fact, these functions are related by the fundamental equation




A related function is the Weierstrass function, which is defined by


The function is not periodic and hence not elliptic. It is holomorphic on , however, and has simple poles at the points of . Furthermore, it is related to the Weierstrass function by the rule


This function arises naturally when calculating the integral in equation (4) for the Weierstrass data of Costa’s surface. Although the function is not periodic, it does satisfy the following period relations:


Furthermore, a theorem due to Legendre states that


These facts are used in the next subsection.

We refer to special points in or as marked points. These are points at which singularities, such as poles, occur. Open disks centered at marked points are referred to as marked disks. Let denote the disk of radius with center at the point , and define by

The set can be viewed as the collection of all points in that are mapped to a marked disk in by the projection .

We now specialize to the case where and . In this case, is the standard square lattice of Gaussian integers in , and we simply write for and for . The numbers in equation (10) now satisfy


The set defines a fundamental domain of the covering . Let

where , , and are small positive numbers, and define to be . Note that is a unit square in with four quarter-disks of radius removed from the corners of , two half-disks of radius removed from the midpoints of the horizontal edges of , and two half-disks of radius removed from the midpoints of the vertical edges of . The projection is a torus with three marked disks removed. These disks are centered at , , and .

With these preliminaries, we are ready to specify the Weierstrass data for Costa’s surface. We take our domain to be , and set


Both of these functions are holomorphic on , and since is simply connected, the integration in equation (4) is path independent.

Costa’s Minimal Surface

In general, one might have to resort to numerical integration in equation (4) to obtain the Weierstrass representation for a surface. However, the integration can be carried out explicitly for the functions in equation (16), when restricted to . This calculation was first performed by Alfred Gray and is given in [2]. The result is as follows:

Theorem. Let denote the lattice of Gaussian half-integers. Then the Weierstrass data in equation (16), when substituted into equation (3) and integrated, yields the holomorphic curve , whose components are given by

The corresponding trace of is Costas minimal surface.

We now use the properties of and discussed earlier to demonstrate that Costa’s surface is topologically a torus with three points removed. Define three points in by , , and . We may think of these points as the projection to of in . From the form of in equation (17), it is clear that is holomorphic on .

Note that, for the basic periods , we have as sets. In this case, it is clear from the definition of the function in equation (11) that

This fact, together with equation (14), implies that


Equations (13, 17, 18) then allow us to conclude that

This clearly implies that for all , that is, that is periodic. Thus, we have demonstrated the following:

Proposition: The map , where is given by equation (17), solves the period problem. Hence projects to a real-analytic map .

Substituting equation (16) into equation (6), we see, with the help of equations (9, 15), that the metric on is given by


From this formula, it is not difficult to show that the function has poles of order two at , , and . Thus, the metric diverges at the ends of Costa’s surface, which we define to be the image under of the punctured disks , , and in . Here, denotes the disk in , excluding its center. We shall see later that diverges at different rates at than at and . This is due to the fact that the principal part of the Laurent expansion of at has a larger leading coefficient than that at or , even though all of the poles are of the same order. As can be seen in Figure 3, the end corresponding to the punctured disk about is asymptotically planar, while the other two ends are asymptotically catenoidal in form.

In the next section, we describe the contents of the Minimal Surfaces package. After that, we show how to use the package to generate polyhedral data representing the trace of the parametrization


where Omega^^(epsilon_1,epsilon_2,epsilon_3):=pi[Omega(epsilon_1,epsilon_2,epsilon_3)].

The Minimal Surfaces Package

We now describe the public functions in the Minimal Surfaces package. Note that there are a number of utility functions with private context that are not described here. Further documentation may be found within the package source file.

Special Data Types

is a two-dimensional graphics primitive specifying a circular arc with center , radius , and initial and terminal angles , which are assumed to satisfy . The orientation may be either Clockwise or CounterClockwise.

Surface-Generating Functions

returns the image in corresponding to the point in with coordinates .

returns the normal displacement at distance from along the normal of at , where is a regular map on some domain into .

Vertex Creation and Triangulation

produces a collection of planar points within the boundary specified by the closed curve , which consists of a list of line segments and circular arcs and is assumed to be traversed counterclockwise.

The output is a list of the form . The interior list consists of planar points properly inside the boundary specified by , while is a list of sublists, each member of which is a list of planar points lying on the corresponding line segment or arc specified by .

Currently, the only supported option is , which specifies the number of sample points to use in the and directions. The defaults are .

glues together a pair of polyhedral surfaces specified by the lists and . These lists are assumed to have the same form as output from the Triangulate function. It is also assumed that the boundary components of and have the same shape.

The output has the form , where is the join of the vertex sublists of and , and is a list of sublists of the form .

produces a triangulation of a planar region. The list is of the same form as the output of CreateVertexData, and , which consists of a list of line segments and circular arcs, defines the boundary of the region to be triangulated.

The output is a list of the form . The elements of are planar points of the form . The elements of are ordered triples of positive integers . The elements of such triples refer to positions in the list of vertices and thereby define triangles. The list is of the form , where each sublist contains positive integers pointing to the position in the list of points that lie on the connected component of the boundary of the triangulated region.

This function has one option, which is of the form . The default value for is , in which case no identifications occur. The other possibility is , where each element of specifies the positions of a pair of edges in that are to be identified, and each element of specifies the positions of a cycle of circular arcs whose endpoints are to be identified.

Remark: The result after identification is assumed to be a Riemann surface, possibly with several marked disks removed.

Graphics Functions

returns a list of 3D graphics directives and primitives describing a polyhedron. The sublist consists of points that are the vertices of a polyhedron, and is a list of sublists of the form , where has the form . Each is a list of graphics directives to be applied to the face list .

Remark: This function could easily be modified to allow the application of graphics directives to edges and vertices.

Import/Export Functions

writes geometric data contained in the list to an ASCII file named . The list is assumed to be in the same form as input data to CreatePolyhedron, except that Mathematica graphics directives are replaced with equivalent directives that are compliant with the export format. The form of the output is specified by . Currently the only supported format is PLY, which stands for the polygon file format

Remark: A desirable enhancement to this function would be to include a POVRAY export format. This would produce object data appropriate for creating scenes with the Persistence of Vision Raytracer program.

The Polygon File Format

We now give a brief summary of the PLY file format, limiting our discussion to those aspects used in the current version of the ExportGraphics function. Further details may be found in [5].

A PLY file has three main parts: a header, a list of vertices, and a list of faces. Listing 1 is an example of a PLY file specifying a cube with faces of various colors. The header contains comments and declarations of data types and their properties. Each line of a comment starts with the token , and data types are specified with the token . The elements necessary for our purposes are and . The declaration of an element must include the number of occurrences of elements of that type in the file. The element’s properties are declared immediately after its declaration. For example, the group of statements

element vertex 8
property float32 x
property float32 y
property float32 z

specifies that the file in question contains data for eight vertices and that each vertex has three properties , which are floating-point numbers representing the coordinates of the vertex.

Each line of the vertex list specifies the coordinates of a vertex. A list of data specifying faces follows the vertex data. Each line of this list gives, in order, the number of vertices in a face and the position in the vertex list corresponding to a vertex of the face. Color attributes of the face can then be given in the form of RGB-color intensities, which are specified by integers from 0 to 255. The orientation of a face can be determined from the order in which its vertices are presented. Note that indexing of the vertices begins at 0. Thus, when 0 occurs in the vertex list of a face in Listing 1, for example, it points to the first vertex, which has coordinates .

format ASCII 1.0
comment a cube                          a comment
element vertex 8                        8 vertices in file
property float32 x                      x coordinate of vertex
property float32 y                      y coordinate of vertex
property float32 z                      z coordinate of vertex
element face 6                          6 faces in file
property list uchar int vertex_indices    vertex incidence list
property uchar red
property uchar green
property uchar blue
end header
0 0 0                                   start of vertex list
0 0 1
0 1 1
0 1 0
1 0 0
1 0 1
1 1 1
1 1 0
4 0 1 2 3 128   0 128                    start of face list
4 7 6 5 4 192 128   0
4 0 4 5 1 128 192   0
4 1 5 6 2 256   0   0
4 2 6 7 3   0 256   0
4 3 7 4 0   0   0 256

Listing 1. A simple PLY file specifying a cube with faces of various colors. The comments in italics to the right are not part of the file.

A Note on Performance

Of the functions listed above, the most time-consuming is Triangulate, which internally calls DelaunayTriangulation. It seemed convenient to use the latter because it is part of the ComputationalGeometry package, which is a standard Mathematica add-on package. Another attractive feature of Delaunay triangulation is that, by design, it produces the most regular triangulation of a planar point set. The most efficient algorithms for Delaunay triangulation have a time complexity of .

The model shown in Figure 4 was created by evaluating Triangulate on a vertex set containing about 60,000 points. This took more than 14 hours on a workstation with dual Opteron 244 processors. (Mathematica Version 5.2 only used one CPU, and the availability of RAM was not an issue.)

There are at least two ways the overall performance of the Minimal Surfaces package could be improved, both of which involve reducing the vertex data used as input to Triangulate. This is discussed further in the subsection on enhancements.

We note that an alternative approach might be to replace the current version of Triangulate with an “adaptive cell division algorithm.” In this scenario, one would start with a sparse vertex set in the given domain , which would be triangulated by some method. Then one would use some function, such as the metric in equation (6), to decide if the triangulation needs to be refined in a neighborhood of any given vertex. The best method of refinement is likely to be midpoint subdivision of any triangle incident to the vertex in question. This may prove to be faster than using Delaunay triangulation.

Generating the Model Data for Costa’s Surface

We now show how the functions described above are used to create a model of Costa’s surface. The essential steps are as follows:

  1. Use the Line and Arc data types to define a curve bounding a suitable domain .
  2. Call CreateVertexData to obtain a collection of points .
  3. Apply Triangulate to in order to obtain its Delaunay triangulation; provide boundary identifications, if appropriate.
  4. Apply a composition of ParallelSurface and CostaSurface to the data generated in step 2 to obtain vertex data for two polyhedral surfaces in that are normal displacements of points lying on Costa’s surface. Note that the faces of the resulting polyhedral surfaces are defined by the same incidence relations resulting from step 3.
  5. Use GlueComponents to create a single volume-bounding polyhedron from the two normally displaced polyhedral surfaces created in step 4.
  6. Call ExportGraphics to produce the PLY file required for printing a solid model.

First, we load the Minimal Surfaces package and define some useful functions. It is available from www.mathematica-journal.com/data/uploads/2010/12/MinimalSurfaces.m. We assume the package is contained in the same directory as this notebook.

The following function defines the curve bounding the region described in the paragraph after equation (15). Recall that the input parameters represent the radii of the three marked disks with centers at , , and , where , as before, denotes the projection . Note that this curve consists of 16 segments and is traversed counterclockwise, hence the circular arcs it contains are traversed clockwise.

The function grid generates a rectangular array of horizontal and vertical lines with offset , step size , and index range ; it is used to generate the background grid in Figure 2.

Example Output from CreateVertexData and Triangulate

Before invoking the commands that produce the graphics object illustrated in Figure 3, we illustrate what the output of the CreateVertexData and Triangulate functions looks like for a small mesh size.

The list v produced is of length 2. The first part of v consists of the interior points shown in green in Figure 1. Here is a short listing.

The second part of v consists of 16 sublists, each of which contains vertices that lie on the corresponding segment of b. For example, the last segment in b corresponds to the quarter-circle in the lower-left corner of Figure 1. The vertices that lie in the segment are the elements of the list .

Here is the first part of , which contains vertices that lie on the leftmost horizontal line segment on the real axis. Note that the endpoints of incident segments overlap.

Figure 1, below, graphically illustrates the structure of v. The green points are vertices that lie in the interior of the boundary curve b, and the red points are vertices lying in b itself. Note that all of the points in v lie in the gray rectangular grid.

Figure 1. Vertex data generated using CreateVertexData with parameter values , , and .

We now triangulate the vertex data v shown in Figure 1 without edge identifications. The resulting triangulation differs from a Delaunay triangulation in that certain triangles incident to boundary points of but not lying within the domain itself are excluded.

The output p of Triangulate has three parts. The first part is a flat list of all the vertices (both interior and boundary) shown in Figure 1. Duplicate points have been removed.

The second part of p contains incidence relations that define faces, which are always triangles. The positive integers in each sublist of are locations of vertices in .

For example, the following command extracts the vertices in that correspond to the first face listed in .

The third part of p is a list containing one sublist that defines the boundary curve of the triangulation. When identifications are used, the list may contain multiple boundary components. In the next subsection, for example, Triangulate is called with identifications that produce three boundary components. Each component in that case corresponds to the boundary of a marked disk that has been removed from the torus described in the paragraph after equation (17).

The following group of commands first maps the vertices of p into the horizontal plane of ; the resulting output list is named q. Then some graphics directives are added to q, and the resulting list is used as input to CreatePolyhedron. The output list s is then flattened so that it can be used as input to Graphics3D.

In general, the length of the output s of CreatePolyhedron is the same as the length of the second part of the input qw. Thus, in the present case, s has length 1. The first part of consists of graphics directives, and the second part consists of graphics primitives that are applied to the faces of the triangulation q. Part of the output from the current evaluation is shown below. To create the polyhedron shown in Figure 3, we will build a list similar to qw, above. The second part of the resulting list itself has five parts, each of which corresponds to one of the five parts of the surface shown in the figure, namely, the green (or inside) part, the yellow (or outside) part, and the three purple rims.

Figure 2 shows what the triangulation of the small dataset looks like.

Figure 2. Triangulation of the vertex data in Figure 1 using the Triangulate function.

Creation and Triangulation of Vertex Data for a Torus

In this subsection, we discuss how the CreateVertexData and Triangulate functions were used to help produce the polyhedron shown in Figure 3. The mesh size of 50×50 employed here is moderate in order to reduce computation time and the size of the resulting datasets. Note that our marked disks are given radii of and . This compensates for the difference in the rate of divergence at the poles , , and of the parametrization defined in the proposition following equation (17). The ends of the resulting surface are roughly the same size.

The list edgeIds specifies edges in b that are to be identified, and the list vertexIds specifies the endpoints of circular arcs in b that are to be identified. For example, the sublist in edgeIds indicates that the lower-left horizontal line segment in Figure 2 is to be identified with the upper-left horizontal line segment. Similarly, the sublist in vertexIds indicates that corresponding endpoints of the semicircular arcs with centers at and are to be identified. The list vertexIds is used to identify which components of the boundary curve close up to form boundary curves. In this case, we get three such curves.

When we pass these identifications to Triangulate as an option, output similar to that in the previous subsection is produced, except that some incidence relations in p are reassigned, and extra vertices are dropped.

Vertex Data for Normal Displacements of Costa’s Surface

The ParallelSurface function is used to create a small normal offset of the vector-valued function it takes as an argument. Here, the function is CostaSurface, which can be used to generate the true Costa surface. Note that, since CostaSurface returns a vector in , ParallelSurface also returns a vector in .

To create the polyhedron in Figure 3, we choose as our offset and define two real vector-valued functions accordingly.

The following commands use and to map the vertex data of the two-dimensional triangulation defined by p into .

The function GlueComponents is now used to connect the boundary components of q1 and q2. First, the vertex and face lists of q1 and q2 are concatenated, and then additional faces, which are incident to vertices in the boundary curves, are appended to the face list. The resulting polyhedron bounds a volume in .

We now define the graphics directives we wish to apply to each component of our surface.

The following command produces input in the form required by CreatePolyhedron. That is, qw has the form , where each is a list of incidence relations in the vertex list defining faces and each is a list of graphics directives to be applied to .

Now we apply CreatePolyhedron to get a graphics object. As before, we need to flatten the list s to produce input acceptable to the Graphics3D function.

Figure 3. Costa’s surface with , , , and .

Printing Technology and the Solid Model

The photograph shown in Figure 4 is a 3D model of Costa’s surface that was printed using a ZCorp Model 402Z 3D printer. This device is no longer in production but is similar to the ZCorp Model 510. The principal difference between the two is that the Model 510 has a higher print resolution. See [6] for more details on the Model 510.

We first convert our choice of RGB values to integers in the range .

Now we produce a list rw in the correct format for input into the ExportGraphics function. Note its similarity in structure to qw.

The PLY file corresponding to the surface pictured in Figure 3 is produced by invoking the ExportGraphics function as follows.

To obtain the solid model illustrated in Figure 4, a PLY file was created using the parameters , , , and . The large mesh size ensured that the model would be smooth, and the rather large normal displacement ensured that it would be thick enough to avoid breakage during production. The resulting object is about eight inches in diameter.

Figure 4. Photograph of a solid model of Costa’s surface generated from a PLY file with parameters , , , and .


In this section, we put what was done here in perspective and outline some ways of generating Weierstrass data that may lead to new examples of minimal surfaces. These methods, together with other enhancements discussed in this section, may be incorporated in future versions of the Minimal Surfaces package.

Genus One Minimal Surfaces

First, it would be interesting to look for other examples of genus one minimal surfaces, possibly with more than three ends. The Weierstrass data for such a surface would be given by elliptic functions on . An obvious choice for Weierstrass data would be


where is a constant, and , for , is defined on by the absolutely convergent series

The case is just that of Costa’s surface, where and . The convergence of this series is actually rather subtle for , and the difference in the summands cannot be separated. For , we may rewrite as

where is the Eisenstein series of weight :

It was mentioned earlier that any elliptic function may be expressed as a rational function of and its derivative . We can achieve this for the Weierstrass data in equation (21) by repeatedly differentiating equation (9) and applying the above observations.

Higher Genus Minimal Surfaces

We now discuss two ways in which one may generate examples of minimal surfaces that have a topology other than that of a torus.

Automorphic Functions on Hyperbolic Space

Every compact, connected surface has a simply-connected universal covering space, which we denote by . If has no additional structure, we know that must be topologically a sphere, which we take to be the extended complex plane , or an open subset of the plane, which we take to be . A surface is said to have a conformal structure if, for every point , one has a notion of the angle between any two tangent vectors based at . If is a simply connected region in , we have the following fundamental fact from complex analysis:

Riemann Mapping Theorem. Any simply connected region in the complex plane (other than itself) is conformally equivalent to the unit disk .

A Riemann surface is a surface endowed with a conformal structure. If is a Riemann surface of genus at least two, the Riemann mapping theorem allows us to think of its universal covering space as being the unit disk . In fact, following [7], a simply connected Riemann surface is classified as being

  • elliptic if it is conformally equivalent to the whole Riemann sphere
  • parabolic if it is conformally equivalent to the finite plane
  • hyperbolic if it is conformally equivalent to the unit disk

Thus, any compact Riemann surface is of the form , where is a discrete subgroup of a group of conformal transformations on .

  • If is elliptic, then is the full group of linear fractional transformations for complex constants , , , and with .
  • If is parabolic, then is the group of linear transformations , for complex constants and .
  • If is hyperbolic, then is the group of linear fractional transformations of the form , where , , and is the complex conjugate of . In this case, is isomorphic to the matrix group .

In accordance with the above, we say that is elliptic if is elliptic, is parabolic if is parabolic, etc. To simplify the discussion a bit, we sometimes use the term Weierstrass data to mean a pair of meromorphic functions on a compact Riemann surface . These functions lift naturally to a pair of -invariant functions on .

The catenoid is an example of a minimal surface that arises from Weierstrass data on an elliptic Riemann surface. The main subject of this paper has been Costa’s surface, which arises from Weierstrass data on a parabolic Riemann surface (i.e., a torus). Minimal surfaces of higher topological type arise from Weierstrass data on Riemann surfaces of hyperbolic type. In order to find such data, we seek meromorphic functions on the unit disk that are invariant under a discrete subgroup
. Such functions are referred to as automorphic functions. A fundamental domain of acting on can be specified by a subregion of bounded by line segments and circular arcs. Thus, the rendering methodology used here should extend directly to the hyperbolic case. This leads us to pose the following:

Problem: Use automorphic functions on the unit disk as Weierstrass data for the construction of new examples of minimal surfaces.

It should be noted that the author has not been able to find any references in the mathematical literature taking this approach to the construction of minimal surfaces. There is another approach, however, which is briefly discussed in the next section.

Algebraic Curves in

Riemann surfaces may also be realized as affine algebraic curves in , which are the solutions of polynomial equations in two complex variables of the form . Such curves are not compact, but they can be compactified via an embedding into the complex projective plane . From this perspective, Riemann surfaces can easily be identified with branched coverings of the complex plane or the Riemann sphere . In equation (9), for example, we see that the Weierstrass data for Costa’s surface essentially provides a complex parametrization of the elliptic curve

where and .

In [8], Thayer describes a family of higher genus minimal surfaces, which he calls Chen-Gackstatter-Karcher surfaces or CGK surfaces. In constructing these surfaces, he first identifies Riemann surfaces of the form

where is a rational function in of a specific form. He then constructs minimal surfaces with the Weierstrass data and , which are defined on a specific branch of the underlying Riemann surface. Here, the Gauss map may be viewed as a multivalued function from into the Riemann sphere . This leads us to pose the following:

Problem: Is it possible to construct an explicit correspondence between automorphic functions on the unit disk and affine or projective algebraic curves? In particular, can we find automorphic functions on that are Weierstrass data for CGK surfaces?

We note that Thayer’s paper [8] contains numerous figures illustrating specific examples of CGK surfaces. The figures were generated using a program called MESH, which was written by James Hoffman. MESH is an adaptive mesh generation program that is specifically designed to generate vertex, edge, and face data for minimal surfaces. This data is generated via numerical integration of the component functions of the Weierstrass representation. MESH is a command-line application with a client-server architecture, the use of which requires some knowledge of C++ or FORTRAN programming. See [9] for further documentation and program download.

Periodic Minimal Surfaces

Finally, we note that there are known examples of surfaces that are periodic in the sense that they are invariant under a discrete group of rigid motions in . A particularly interesting family of such surfaces is described in [10]. It would be interesting to search for new examples of such periodic minimal surfaces via the use of automorphic functions.

Enhancements to the Minimal Surfaces Package

We finish with a short list of possible enhancements to the Minimal Surfaces package.

  • Improve performance by employing a method for selecting the density of vertex data on the basis of the surface metric given by equation (6). Alternatively, use the metric as a basis for an adaptive mesh generation algorithm.
  • Improve performance by employing symmetries of surfaces to reduce the amount of computation.
  • Include functionality for rendering other minimal surfaces, including the CGK surfaces described above.
  • Incorporate functionality that permits the visualization of different coordinate systems on a minimal surface. Of particular interest would be principal curvature lines and asymptotic lines.
  • Include functionality for rendering geodesics, including the ability to find and render closed geodesics.
  • Include an export function that is suitable for use with the Persistence of Vision Ray Tracer (POVRAY) program.

Finally, we note that the Minimal Surfaces package was originally written using Mathematica Version 5.2. One reviewer of this paper indicated that the DelaunayTriangulation routine in the Computational Geometry Package for Mathematica 5.2 is known to have poor complexity and has suggested the use of Martin Kraus’s Polygon Triangulation package [11] instead. Martin Kraus’s package seems to provide significantly better performance and may be incorporated into a future version of Minimal Surfaces. The same reviewer also pointed out that, in later releases of Mathematica, the built-in Export function supports both PLY and POVRAY file formats.


The idea for this paper grew out of a short course on computational topology the author attended from July 6-16, 2004. The course was taught by Herbert Edelsbruner (Duke University) and John Harer (Duke University), and held at the Institute for Mathematics and Its Applications (IMA) at the University of Minnesota. Thanks are due to Doug Arnold, the director of the IMA, for hosting the course, Herbert Edelsbruner for kindly providing access to his 3D printer, and Rachael Brady (Duke University) for her patient help in debugging the export function of the Minimal Surfaces package and successfully printing out the solid model shown in Figure 4. Thanks are also due to the reviewers of this paper, who provided many useful comments.


[1] R. S. Palais, “The Visualization of Mathematics: Towards a Mathematical Exploratorium,” Notices of the American Mathematical Society, 46(6), 1999 pp. 647-658.
[2] A. Gray, Modern Differential Geometry of Curves and Surfaces with Mathematica, 2nd ed., Boca Raton, FL: CRC Press, 1999.
[3] D. A. Hoffman and W. Meeks III, “A Complete Embedded Minimal Surface in with Genus One and Three Ends,” Journal of Differential Geometry, 21, 1985 pp. 109-127.
[4] K. Chandrasekharan, Elliptic Functions, Grundlehren der Mathematischen Wissenschaften 281, New York: Springer-Verlag, 1985.
[5] P. Bourke. “PLY-Polygon File Format.” local.wasp.uwa.edu.au/~pbourke/dataformats/ply.
[6] ZCorporation. www.zcorp.com.
[7] G. Springer, Introduction to Riemann Surfaces, 2nd. ed., New York: Chelsea, 1981.
[8] E. C. Thayer, “Higher-Genus Chen-Gackstatter Surfaces and the Weierstrass Representation for Surfaces of Infinite Genus,” Experimental Mathematics, 4(1), 1995 pp. 19-39. www.emis.de/journals/EM/expmath/volumes/4/4.html.
[9] J. T. Hoffman, “MESH: A Program for Generating Parametric Surfaces Using an Adaptive Mesh,” 1996.
[10] V. R. Batista, “A Family of Triply Periodic Costa Surfaces,” Pacific Journal of Mathematics, 212(2), 2003 pp. 347-370. citeseerx.ist.psu.edu/viewdoc/summary?doi=
[11] M. Kraus. “Polygon Triangulation.” library.wolfram.com/infocenter/MathSource/23/.
O. Michael Melko, “Visualizing Minimal Surfaces,” The Mathematica Journal, 2010. dx.doi.org/doi:10.3888/tmj.12-6.

About the Author

Mike Melko is an assistant professor of mathematics at Khalifa University, Sharjah Campus. He received his Ph.D. from the University of California at Santa Cruz in 1989. His research interests are in the areas of differential geometry, mathematical physics, and computational mathematics.

O. Michael Melko
Khalifa University of Science, Technology and Research
Sharjah Campus
PO Box 573