18.3. Mandelbrot Example

C'mon, now, what kind of a book on graphics programming would this be without an example of drawing the Mandelbrot set?

Our last shader of the chapter doesn't fall into the category of attempting to achieve an artistic or painterly effect, but it's an example of performing a type of general-purpose computation on the graphics hardware for scientific visualization. In this case, the graphics hardware enables us to do real-time exploration of a famous mathematical function that computes the Mandelbrot set. In 1998, Michael Rivero created a RenderMan shader that creates the Mandelbrot set as a texture. Subsequently, Dave Baldwin adapted this shader for GLSL, Steve Koren then made some modifications to it to get it working on programmable hardware for the first time, and I've adapted this shader further for inclusion in this book.

The point of including this shader is to emphasize the fact that with the OpenGL Shading Language, computations that were previously possible only in the realm of the CPU can now be executed on the graphics hardware. Performing the computations completely on the graphics hardware is a big performance win because the computations can be carried out on many pixels in parallel. Visualization of a complex mathematical function such as the Mandelbrot set barely begins to tap the potential for using programmable shader technology for scientific visualization.

18.3.1. About the Mandelbrot Set

To understand the Mandelbrot set, we must first recall a few things from our high school mathematics days. No real number, when multiplied by itself, yields a negative value. But the imaginary number called i is defined to be equal to the square root of 1. With i, the square root of any negative number can easily be described. For instance, 3i squared is 9, and conversely the square root of 9 is 3i.

Numbers that consist of a real number and an imaginary number, such as 6 + 4i, are called complex numbers. Arithmetic operations can be performed on complex numbers just as on real numbers. The result of multiplying two complex numbers is as follows:

x = a + bi
y = c + di
xy = ac + adi + cbi - bd
= (ac - bd) + (ad + bc)i

Because complex numbers contain two parts, the set of complex numbers is two-dimensional. Complex numbers can be plotted on the complex number plane, which uses the horizontal axis to represent the real part and the vertical axis to represent the imaginary part (see Figure 18.6).

In Figure 18.6, we see three symbols plotted on the complex number plane. A small square is plotted at the complex number 2 + i, a small circle is plotted at 1 + 2i, and a triangle is plotted at 1.5 0.5i.

With the assistance of 1970s computer technology (quite inferior by today's consumer PC standards), a mathematician named Benoit Mandelbrot began studying a recursive function involving complex numbers:

In this function, the value of Z is initialized to c, the value of the complex number being tested. For each iteration, the value of Z is squared and added to c to determine a new value for Z. This amazingly simple iterative formula produces what has been called the most complex object in mathematics, and perhaps the most complex structure ever seen!

It turns out that for some values of c, the function eventually approaches infinity, and for other values it does not. Quite simply, values of c that go off to infinity are not part of the Mandelbrot set, and the rest are. If we use a computer (or an OpenGL Shading Language-capable graphics accelerator!) to test thousands of points on the complex number plane and assign a gray level to those that go off to infinity and black to those that don't, we see the following familiar cardioid/prickly pear shape start to appear (see Figure 18.7).

How exactly do we test whether these values go off to infinity? Well, Mandelbrot helped us out a bit here. He showed that if the magnitude of Z (its distance from the origin) was greater than 2, the value of the function would go to infinity. To encode this function in a programming language, all we need to do is stop iterating when the magnitude of Z surpasses 2. Even easier, because we're always dealing with Z², we can simply check to see if Z² is greater than 4.

The values inside the black region of Figure 18.7 do not go off to infinity in any reasonable number of iterations. How do we deal with these? To prevent our computer or graphics hardware from locking up in an infinite loop, we need to decide on the maximum number of iterations to allow before we give up and assume the point is inside the Mandelbrot set. With these two exit criteria, we can write a simple loop that computes the Mandelbrot set.

The beauty of the Mandelbrot set can be enhanced if we color code the number of iterations needed to determine whether a particular point is inside the set. Values determined to be outside the set on the first iteration are given one color value, values determined to be outside the set on the second iteration are given another color value, and so on. In Figure 18.7, I've used gray levels to indicate the number of iterations. The medium gray on the outside represents values that are identified as outside the Mandelbrot set on the very first iteration. Values in white along the edge took 20 iterations to be identified as being outside the Mandelbrot set.

The edges of the Mandelbrot set hold an infinite amount of self-similar variety. By continuing to zoom in on the edge detail, you can find complex numbers whose magnitude stayed below 2 for hundreds or even thousands of iterations of the function before finally exceeding the threshold and zooming off to infinity.

Here are a few other deep thoughts that you can use to amaze and amuse your friends:

• The length of the border for the Mandelbrot set is infinite.

• All the regions inside the Mandelbrot set (i.e., the black regions) are connected.

• There is exactly one band surrounding the Mandelbrot set for each iteration value (e.g., a band that exceeded the threshold on the first iteration, a band that exceeded on the second iteration, and so on). The iteration bands go completely around the Mandelbrot set, do not break, and do not cross each other. When you've zoomed in to explore an edge region with amazing complexity, this is pretty astonishing.

• There are an infinite number of "mini-Mandelbrots" (regions that look like warped or transformed versions of the full Mandelbrot set) within the original.

18.3.2. Vertex Shader

The vertex shader for the Mandelbrot set (see Listing 18.6) is almost exactly the same as the vertex shader for the simple brick example that was described in Section 6.2. The only difference is that we assume that texture coordinates will be provided in the range of [0,1] for both s and t, and we map these values into the range [-2.5, 2.5] and store the result into a varying variable named Position. This enables the fragment shader to plot values directly onto a coordinate system that is just the right size for plotting the Mandelbrot set, and it has the point (0,0) in the middle. If the application draws a single polygon that is a screen-aligned square and has texture coordinate (0,0) in the lower-left corner and (1,1) in the upper right, the result is a standard representation of the Mandelbrot set. Of course, with our OpenGL Mandelbrot shader, the Mandelbrot set can be "textured" onto any geometry and we even apply a simple lighting model as an added bonus.

Listing 18.6. Vertex shader for drawing the Mandelbrot set

uniform vec3 LightPosition; uniform float SpecularContribution; uniform float DiffuseContribution; uniform float Shininess; varying float LightIntensity; varying vec3 Position; void main() { vec3 ecPosition = vec3(gl_ModelViewMatrix * gl_Vertex); vec3 tnorm = normalize(gl_NormalMatrix * gl_Normal); vec3 lightVec = normalize(LightPosition - ecPosition); vec3 reflectVec = reflect(-lightVec, tnorm); vec3 viewVec = normalize(-ecPosition); float spec = max(dot(reflectVec, viewVec), 0.0); spec = pow(spec, Shininess); LightIntensity = DiffuseContribution * max(dot(lightVec, tnorm), 0.0) + SpecularContribution * spec; Position = vec3(gl_MultiTexCoord0 - 0.5) * 5.0; gl_Position = ftransform(); }

18.3.3. Fragment Shader

The fragment shader implements the algorithm described in the previous section. Uniform variables establish the maximum number of iterations and the starting point for viewing the Mandelbrot set (center value and zoom factor). The application is given artistic license to use uniform variables to set one color for points inside the set and two colors to use for points outside the set. For values outside the set, the color gradient from OuterColor1 to OuterColor2 is broken into 20 separate bands, and the cycle is repeated if the number of iterations goes above 20. It is repeated again if the number of iterations goes above 40, and so on.

This shader maps the x coordinate of the computed position in the complex number plane (i.e., the value in Position.x) to the real number in the iterative function, and the y coordinate to the imaginary number. After the initial conditions have been established, the shader enters a loop with two exit criteriaif we reach the maximum number of iterations allowed or if the point proves to be outside the set. Within the loop, the function Z² + c is computed for use in the next iteration. After the loop is exited, we compute the color of the fragment. If we're inside the set, we use the inside color. If we're on the edge or outside, we blend the edge color and the outer color, depending on the number of iterations that have occurred.

The complete fragment shader is shown in Listing 18.7.

Listing 18.7. Fragment shader for computing the Mandelbrot set

varying vec3 Position; varying float LightIntensity; uniform float MaxIterations; uniform float Zoom; uniform float Xcenter; uniform float Ycenter; uniform vec3 InnerColor; uniform vec3 OuterColor1; uniform vec3 OuterColor2; void main() { float real = Position.x * Zoom + Xcenter; float imag = Position.y * Zoom + Ycenter; float Creal = real; // Change this line. . . float Cimag = imag; // . . .and this one to get a Julia set float r2 = 0.0; float iter; for (iter = 0.0; iter < MaxIterations && r2 < 4.0; ++iter) { float tempreal = real; real = (tempreal * tempreal) - (imag * imag) + Creal; imag = 2.0 * tempreal * imag + Cimag; r2 = (real * real) + (imag * imag); } // Base the color on the number of iterations vec3 color; if (r2 < 4.0) color = InnerColor; else color = mix(OuterColor1, OuterColor2, fract(iter * 0.05)); color *= LightIntensity; gl_FragColor = vec4(color, 1.0); }

There is obviously room for improvement in this shader. One thing you might do is improve the color selection algorithm. One possibility is to use a 1D texture to store a color lookup table. The number of iterations could be used to index into this table to obtain the color for drawing the fragment.

After you've invented a pleasing coloring scheme, you can explore some of the popular Mandelbrot "tourist locations." Various books and Web sites have published the coordinates of interesting locations in the Mandelbrot set, and these shaders are set up so that you can plug those coordinates in directly and zoom in and see for yourself. Figure 18.8 shows a few that I explored.

18.3.4. Julia Sets

Julia sets are related to the Mandelbrot set. Each point in the Mandelbrot set can be used to generate a Julia set, and these are just as much fun to explore as the Mandelbrot set. The only difference is that the constant c in the equation term Z² + c is initialized to the value of a point in the Mandelbrot set other than the one currently being plotted. To change the Mandelbrot shader into a fragment shader for doing Julia sets, change the two lines of code that initialize the value of c:

float  Creal = -1.36;     // Now we'll see an interesting Julia set
float  Cimag = 0.11;

Figure 18.9 shows a few examples of the Julia sets that can be rendered with this shader. You might want to change this shader so that the values for Creal and Cimag can be passed in as uniform variables when you want to draw a Julia set. The numbers are no longer imaginary. Now we can do real-time exploration of the mathematical universe through the OpenGL Shading Language!