Who says math is not beautiful? Anyone who doubts the beauty in math must check out algebraic surfaces.

An implicit functions has the form:

where,  is an arbitrary constant.

For example, the implicit equation of the unit circle is  and the equation  describes a sphere of radius 1.

An algebraic surface is described by an implicit function . It may be rendered using Mayavi by evaluating the expression  on a regular (specified) grid, and then generate isosurface of contour value equal to 0.

Here are some examples of algebraic surfaces plotted using Mayavi:

I wrote a quick function called implicit_plot for plotting algebraic surfaces using Mayavi. The most important argument to the function is of course the expression   as a string. It is probably not the best function, but the idea is to show how to plot implicit surfaces. The code snippet is included at the end of this post. Suggestions for improvement are always welcome.

Let’s start with a simple sphere in three dimensional space, whose equation is  . We can call the function implicit_plot like so (please note that it is assumed that mayavi has been imported in the script as mlab ):

import mayavi.mlab as mlab
figw = mlab.figure(1, bgcolor=(0.1, 0.1, 0.1), size=(400, 400))
implicit_plot('x**2 + y**2 + z**2 - {R:f}**2'.format(R=1), (-3, 3, -3, 3, -3, 3),
              fig_handle=figw, Nx=64, Ny=64, Nz=64, col_isurf=(0.0,0.2,0.8),
              opaque=True, ori_axis=False)
mlab.show()

Running the above code should render a sphere of unit radius in the given Mayavi figure window:

The equation represents a sphere when  and  is the square of the radius.

Here we study the various forms of surfaces represented by the above equation for ;  (i.e. for even )

figw = mlab.figure(1, bgcolor=(0.97, 0.97, 0.97), size=(400, 400))
exponent = [2, 4, 6, 8, 10]
num_exponent = len(exponent)
col_isurf_arr = [(0.0,0.2,0.8), (0.2,0.5,0.4), (0.4,0.8,0.2),
                 (0.8,0.7,0.1), (1.0,0.2,0.0)] # colors
b = 100.0
for i, ex in enumerate(exponent):
    implicit_plot('x**{n} + y**{n} + z**{n} - {b}'.format(n=ex, b=b),
                  (-15, 15, -15, 15, -15, 15), fig_handle=figw,
                  Nx=64, Ny=64, Nz=64, col_isurf=col_isurf_arr[i],
                  opa_val=0.35+0.25*i/num_exponent, opaque=False, ori_axis=True)

mlab.show()

which produces the following output:

We can see that when is even, the sphere transforms into a cube-like surface as the value of increases. The following surface is generated for :

figw = mlab.figure(1, bgcolor=(0.1, 0.1, 0.1), size=(400, 400))
implicit_plot('x**20 + y**20 + z**20 - 100',
              (-2, 2, -2, 2, -2, 2), col_osurf=(0.87,0.086,0.086),
              fig_handle=figw, opa_val=1.0, opaque=True, ori_axis=False)
mlab.show()

Next, we can study the surface transformations for odd values of . i.e. for for . The figure also renders a sphere (exponent=2) for reference.

figw = mlab.figure(1, bgcolor=(0.1, 0.1, 0.1), size=(400, 400))
exponent = [2, 3, 5, 7, 9]
num_exponent = len(exponent)
col_isurf_arr = [(0.0,0.2,0.8), (0.2,0.5,0.4), (0.4,0.8,0.2),
                 (0.8,0.7,0.1), (1.0,0.2,0.0)] # colors
b = 100.0
for i, ex in enumerate(exponent):
    implicit_plot('x**{n} + y**{n} + z**{n} - {b}'.format(n=ex, b=b),
                  (-15, 15, -15, 15, -15, 15), fig_handle=figw,
                  Nx=64, Ny=64, Nz=64, col_isurf=col_isurf_arr[i],
                  opa_val=0.35+0.25*i/num_exponent, opaque=False, ori_axis=True)
mlab.show()

Multiple algebraic surfaces can be rendered by evaluating a product of the expressions. For example, the expression for spheres is .

In the following example we render just two spheres:

figw = mlab.figure(1, bgcolor=(0.1, 0.1, 0.1), size=(400, 400))
# draw a base-plane
implicit_plot('(z + 0)', (-5, 5, -5, 5, -5, 5), fig_handle=figw,
             col_isurf=(0.1,0.1,0.1), opa_val=1.0, opaque=False, ori_axis=True)

# The two spheres
implicit_plot('(x**2 + y**2 + (z-1.414)**2 - 2)*((x-1.5)**2 + (y-1.5)**2 + (z-0.707)**2 - 0.5)',
              (-5, 5, -5, 5, -5, 5), fig_handle=figw, Nx=100, Ny=100, Nz=100,
              col_isurf=(0.87,0.086,0.086), opa_val=1.0, opaque=False, ori_axis=False)
mlab.show()

A 4-way tubing can be constructed by “adding” two cylinder surfaces together. As we saw above, two surfaces may be “added” together by multiplying the two surfaces equations. The equation for the 4-way tubing is thus , where the parameter controls the smoothness of the joint such that increasing increases the smoothness.

figw = mlab.figure(1, bgcolor=(0.1, 0.1, 0.1), size=(400, 400))
implicit_plot('(x**2 + y**2 - 1)*(x**2 + z**2 - 1) - {a}'.format(a=0.005),
              (-5, 5, -5, 5, -5, 5), fig_handle=figw, Nx=201, Ny=201, Nz=201,
              col_isurf=(1.0,204./255,51./255), col_osurf=(1.0,102./255,0.0),
              opaque=True, ori_axis=False)
mlab.show()

Here are a few more examples of some interesting surfaces. The first one is known as the Zitrus surface whose equation is

(since the basic pattern of the code is same, we will just render the figure here)

The next surface is Diabolo, which as the equation of the form

And finally, we render the well known Sweet surface whose expression is

Obviously, the algebraic surfaces are very beautiful, and there are really endless types of them for fun and study. The following extremely resourceful links can help the interested explorers of algebraic surfaces:

1. Implicit function, Wikipedia
2. List of complex and algebraic surfaces
3. Algebraic surfaces Gallery, Hauser Herwig
4. Algebraic surface, Wikipedia
5. Algebraic surfaces, Gerhard Brunthaler
6. Algebraic geometry, Wikipedia
7. Algebraic variety, Wikipedia
8. The Scientific Graphic Project, by David Hoffman, James Hoffman, Matthias Weber, et. al 
9. Surfaces at mathcurve

Here is the code snippet for the implicit_plot function

def implicit_plot(expr, ext_grid, fig_handle=None, Nx=101, Ny=101, Nz=101,
                 col_isurf=(50/255, 199/255, 152/255), col_osurf=(240/255,36/255,87/255),
                 opa_val=0.8, opaque=True, ori_axis=True, **kwargs):
    """Function to plot algebraic surfaces described by implicit equations in Mayavi

    Implicit functions are functions of the form

        `F(x,y,z) = c`

    where `c` is an arbitrary constant.

    Parameters
    ----------
    expr : string
        The expression `F(x,y,z) - c`; e.g. to plot a unit sphere, the `expr` will be `x**2 + y**2 + z**2 - 1`
    ext_grid : 6-tuple
        Tuple denoting the range of `x`, `y` and `z` for grid; it has the form - (xmin, xmax, ymin, ymax, zmin, zmax)
    fig_handle : figure handle (optional)
        If a mayavi figure object is passed, then the surface shall be added to the scene in the given figure. Then, it is the responsibility of the calling function to call mlab.show().
    Nx, Ny, Nz : Integers (optional, preferably odd integers)
        Number of points along each axis. It is recommended to use odd numbers to ensure the calculation of the function at the origin.
    col_isurf : 3-tuple (optional)
        color of inner surface, when double-layered surface is used. This is also the specified color for single-layered surface.
    col_osurf : 3-tuple (optional)
        color of outer surface
    opa_val : float (optional)
        Opacity value (alpha) to use for surface
    opaque : boolean (optional)
        Flag to specify whether the surface should be opaque or not
    ori_axis : boolean
        Flag to specify whether a central axis to draw or not

    """
    if fig_handle==None:  # create a new figure
        fig = mlab.figure(1,bgcolor=(0.97, 0.97, 0.97), fgcolor=(0, 0, 0), size=(800, 800))
    else:
        fig = fig_handle
    xl, xr, yl, yr, zl, zr = ext_grid
    x, y, z = np.mgrid[xl:xr:eval('{}j'.format(Nx)),
                       yl:yr:eval('{}j'.format(Ny)),
                       zl:zr:eval('{}j'.format(Nz))]
    scalars = eval(expr)
    src = mlab.pipeline.scalar_field(x, y, z, scalars)
    if opaque:
        delta = 1.e-5
        opa_val=1.0
    else:
        delta = 0.0
        #col_isurf = col_osurf
    # In order to render different colors to the two sides of the algebraic surface,
    # the function plots two contour3d surfaces at a "distance" of delta from the value
    # of the solution.
    # the second surface (contour3d) is only drawn if the algebraic surface is specified
    # to be opaque.
    cont1 = mlab.pipeline.iso_surface(src, color=col_isurf, contours=[0-delta],
                                      transparent=False, opacity=opa_val)
    cont1.compute_normals = False # for some reasons, setting this to true actually cause
                                  # more unevenness on the surface, instead of more smooth
    if opaque: # the outer surface is specular, the inner surface is not
        cont2 = mlab.pipeline.iso_surface(src, color=col_osurf, contours=[0+delta],
                                          transparent=False, opacity=opa_val)
        cont2.compute_normals = False
        cont1.actor.property.backface_culling = True
        cont2.actor.property.frontface_culling = True
        cont2.actor.property.specular = 0.2 #0.4 #0.8
        cont2.actor.property.specular_power = 55.0 #15.0
    else:  # make the surface (the only surface) specular
        cont1.actor.property.specular = 0.2 #0.4 #0.8
        cont1.actor.property.specular_power = 55.0 #15.0

    # Scene lights (4 lights are used)
    engine = mlab.get_engine()
    scene = engine.current_scene
    cam_light_azimuth = [78, -57, 0, 0]
    cam_light_elevation = [8, 8, 40, -60]
    cam_light_intensity = [0.72, 0.48, 0.60, 0.20]
    for i in range(4):
        camlight = scene.scene.light_manager.lights[i]
        camlight.activate = True
        camlight.azimuth = cam_light_azimuth[i]
        camlight.elevation = cam_light_elevation[i]
        camlight.intensity = cam_light_intensity[i]
    # axis through the origin
    if ori_axis:
        len_caxis = int(1.05*np.max(np.abs(np.array(ext_grid))))
        caxis = mlab.points3d(0.0, 0.0, 0.0, len_caxis, mode='axes',color=(0.15,0.15,0.15),
                              line_width=1.0, scale_factor=1.,opacity=1.0)
        caxis.actor.property.lighting = False
    # if no figure is passed, the function will create a figure.
    if fig_handle==None:
        # Setting camera
        cam = fig.scene.camera
        cam.elevation(-20)
        cam.zoom(1.0) # zoom should always be in the end.
        mlab.show()

Please note that the main block in the the above function is really the following three lines, which generates a scalar field on regular grid using the Mayavi pipleline, and then generates a zero-value isosurface:

scalars = eval(expr)
src = mlab.pipeline.scalar_field(x, y, z, scalars)
cont1 = mlab.pipeline.iso_surface(src, color=col_isurf, contours=[0],transparent=False, opacity=opa_val)

Hope you enjoyed the post.