Inside Koval's 3D Grapher

Koval's 3D Grapher is an interactive function plotter in three dimensions. Here's an example plot:

Koval's 3D Grapher is a conglomerate of various mathematical identities and algorithms and a large chunk of good 'ol hacking. What follows is a summary of some of the main components, their challenges, and how they evolved over time.

Plotting Points

A surface \(f(x,y)\) is made up of an infinite number of points. A sphere is made by drawing every point that is a certain distance from another point. An sinusoidal surface is made up of an infinite number of points that form sine and cosine cross sections. Even a degenerate surface, like a line, is made up of points.

Because a function is a relation between three variables (\(x\), \(y\) and \(z\)), all surfaces are guaranteed to be made up of points that have a distinct value for these variables. Furthermore, these points satisfy the given equation: \[g(x,y,z)=0\] For the unit sphere, the euclidean distance from every point to the origin is one, so \(g(x,y,z)=\sqrt{x^2+y^2+z^2}-1=0\). For a sinusoid, we have \(g(x,y,z)=\sin{x}+\sin{y}-z=0\)

Each point fulfills the equation \(\sin{x}+\sin{y}-z=0\)

If we were to try to find these points, we might naively attempt to iterate through a bunch of \((x,y,z)\) points and see if any of them fulfill our equation. Perhaps we would try points in an orderly manner, and perhaps we would even gauge how "far" a point is from the surface by the magnitude of \(g(x,y,z)\). This works, but it is slow. In order to transverse a region in three dimensions, we must iterate through \(x\), \(y\), and \(z\), which is \(O(n^3)\).

The image on the left shows \(11^3\) points evenly distributed over an interval. The image on the right shows the points for which \(|\sin{x}+\sin{y}-z|<\frac{1}{2}\).

Additionally, actually finding the surface can be difficult, since \(g(x,y,z)\) may get arbitrarily close to zero without actually passing through zero, leading to misplaced points. Consider the relation \( g(x,y,z)=\left( \frac{x}{z} \right) ^2 + \left( \frac{y}{z} \right) ^2 \). Except at \(x=0\) and \(y=0\), \(g\) is never going to equal zero, despite getting arbitrarily close as \(z\to \pm \infty\). How do you decide which points are "close enough" to 0?

This is why, if you've ever worked with 3D functions, you have seen functions in the form \(z=f(x,y)\), \(y=f(x,z)\) or \(x=f(y,z)\). Instead of testing every possible \((x,y,z)\) point (which is of course impossible), we iterate through two variables and solve for the third. Consider our sinusoidal surface. Instead of solving for all our variables in the equation \(\sin{x}+\sin{y}-z=0\), we solve for \(z\) using \(z=\sin{x}+\sin{y}\). We still cannot guarantee that we are going to capture all of the points on the surface (there are an infinite number of them), but we can at least guarantee that all the points will lie exactly on the surface.

The image on the left shows an orthographic overhead view of the surface on the right. Note that the grid corresponds to selected values of \(x\) and \(y\) and the height is computed.

The algorithm to plot points on a surface is relatively simple. Pick an orderly selection of two of the variables, solve for the third, and add the point to a list.

From Points to a Surface

Points are great. They are the basis of our surface. But we can only plot so many points. Eventually, we are going to have to stop making points and start connecting them. To do this, we need to interpolate between points – make our best guess as to where the surface should be. The question is, which interpolation method is best? Which gives the viewer the impression that the entire surface is plotted instead of a selection of points?

How should we interpolate between two points? Is a squiggle better than a straight line?

You've probably guessed that the straight line is the best bet. There are several reasons for this. First of all, lines are quite easy to work with from a graphics perspective. They are computationally cheap and require the least effort to implement. Second, it is much easier to shade a polygon than an arbitrary 2D shape. Since straight lines project onto planes in straight lines, polygons can be colored in 2D on the projected surface, unlike curves. I'll get into this more in Part II. Finally, all surfaces are locally planar, which is the basis of most of 3D calculus. Since we are treating our points as infinitesimal changes in our surface, it makes sense to connect them as if they were actually infinitesimally far away – with straight lines.

Straight line connections make life much easier. However, there still needs to be a way of drawing all these connections in a way that doesn't draw some twice. For reasons that will become apparent in the "Cleaning Up the Mess" section, I chose to use a doubly linked 2D array. In this system, each point is linked to its neighbors, one to the right and left, one above and one below. That way, we can know what connections exist come drawing time.

Each point has as pointer to each of its four neighbors.

It is important to note that the words "above," "below," "right," and "left" have little physical meaning. Especially when we begin to parameterize functions, these connections will have little correlation in the \(xyz\) coordinate world. Think of them as connections to the point plotted before it, the point plotted after it, and the points plotted a row ahead and behind.

This system makes drawing very easy. For every point, draw two lines – one from the point to the point below it and one to the point on the right. This guarantees that no line will be drawn twice, unlike a polygon system.

Each point connected to two of its neighbors can fill an entire surface.

Using this system, we can now plot a rudimentary surface.

A paraboloid surface: \(z=\frac{x^2+y^2}{10}\)

Cleaning Up the Mess

At this point is may seem like the surface is complete. Unfortunately, a surface is often more complex than this – literally. Consider the function \[z=\sqrt{25-x^2-y^2}\] You may recognize this as a half sphere. Perhaps you've done some algebra in your head and noticed that this is nearly equivalent to \(x^2+y^2+z^2=25\). The only difference (and what makes it a half sphere) is that the square root function only returns positive values, so \(z>0\).

See what happens when \(x=5\) and \(y=5\): \[z=\sqrt{25-5^2-5^2}\] \[z=\sqrt{-25}\] \[z=5i\] \(z\) is a complex number. Complex numbers don't fit into the standard Cartesian 3D coordinate system, so \((5,5,5i)\) should not be shown, even though it exists.

This isn't that much of a problem until you realize what happens to the surface after all the complex points are taken out:

Wow! That is truly ugly.

What results is points that are near the edge of the surface – that is, the boundary between real and nonreal points – but not quite there. Because the points are created in a square grid, the points are almost never going to be on a boundary that is any sort of curve, like a circle. This means that there is going to be jagged edges, even with a high number of points.

A solution takes advantage of the fact that the points are linked to their neighbors. Consider a point that is part of the jagged edge seen in the previous picture. One of its neighbors is complex and not shone, but it still exists. Therefore, the real boundary between real and nonreal points exists somewhere between these two points.

A boundary exists somewhere between a real point and its bordering imaginary point.

Thus we need to implement a search algorithm. Instead of an array of values, however, there is an interval of infinite booleans – real or nonreal. Additionally, it is an ordered list since the surface transitions to imaginary then stays there. This is the perfect task for a binary search.

A binary search is also a great time to use recursion. It is essential the reduction of an interval in which the boundary lies, so the process remains the same no matter the size of the interval. In fact, here's the exact code in Koval's 3D Grapher:

function findParametricEdge(realP, nonRealP, xFunc, yFunc, zFunc, iterations){
    //when we have completed the necessary amount of iterations,
    //  return the real end of the interval containing the boundary
    if(iterations === 0){
        return realP;
    }
    
    //find the midpoint of realP and nonRealP in terms of parameterized variables
    var u = (realP.u + nonRealP.u)/2;
    var v = (realP.v + nonRealP.v)/2;
    var scope = {
        u:u,
        v:v,
        x:null,
        y:null,
        z:null
    }
    
    //set any animation variables
    for(var j = 0; j < animationVars.length; j++){
        scope[animationVars[j].name] = animationVars[j].value;
    }
    
    //evaluate function at new point (twice to account for function like x = u, y = x, z = y)
    xFunc.eval(scope);
    yFunc.eval(scope);
    zFunc.eval(scope);
    xFunc.eval(scope);
    yFunc.eval(scope);
    zFunc.eval(scope);
    //create new point containing xyz and uv coordinates for midppoint
    var inBetween = {x:getRealPart(scope.x),y:getRealPart(scope.y),z:getRealPart(scope.z)};
    makeFiniteIfInfinate(inBetween)
    
    //domain stretching for non-cubic domains
    inBetween.x = spreadCoord(inBetween.x, domain.x.min, domain.x.max, domain.x.spread);
    inBetween.y = spreadCoord(inBetween.y, domain.y.min, domain.y.max, domain.y.spread);
    inBetween.z = spreadCoord(inBetween.z, domain.z.min, domain.z.max, domain.z.spread);
    inBetween.u = u;
    inBetween.v = v;
    
    //return new interval
    if(isNonReal(inBetween)){
        return findParametricEdge(realP, inBetween, xFunc, yFunc, zFunc, iterations - 1);
    }else{
        return findParametricEdge(inBetween, nonRealP, xFunc, yFunc, zFunc, iterations - 1);
    }
}    

You'll find out in a later section why everything is defined in terms of \(u\) and \(v\), that is, parameterized. For now we'll stick with the edge detection.

If you run this algorithm on every interval containing both real and imaginary points, then the surface begins to look much more organic – less blocky. Take a look at the same half sphere as a more complete surface:

Looking a lot better.

This is much improved. Now the edge of the surface is clearly seen and most importantly smooth. There is just one thing wrong: loose ends. If you look closely enough you can see that the lines extending to the boundary have only one connection, their base point. This is a large problem when coloring, since coloring depends on polygons to fill in. Without some sort of loop in lines (that is what a polygon is, after all), there is no polygon possible, and so no coloring. What we need is a way to connect all these ends into a cohesive surface boundary, a line that allows polygons to be formed between each end and its neighbor.

To this end, it turns out that simplicity is best. When first mulling over this problem, you might be tempted to try something fancy. One idea is to propagate back from all the loose ends concurrently, traversing edges to extend a loose end "circle" in all directions. When circles from two ends overlap, we connect them.

This is great, but also fiendishly complicated. Also, there are all sorts of edge cases to be worried about. How do you make sure that an end is on the same side of a surface as another end?

In the end, I settled for a shortest distance path algorithm. It proceeds as follows:

1. Make a list of all the ends. This is easily done while plotting.
2. Start with the first end.
3. Find the (unused) end closest in Euclidean distance to this end.
4. Connect these two ends.
5. Repeat steps 3 - 4 until there are no ends left, or the distance between ends is prohibitively far, indicating a new "end loop."

Each dot represents a loose end. The algorithms starts on the left red dot and transverses the left polygon, greedily connecting to the nearest dot. When it gets to the most left dot, the nearest unused dot is farther than the start of the path, so it connects back, creating a loop. The process is then restarted, starting on the right red dot.

And finally, we have our completed surface, ready for coloring.

The final surface, sans coloring, perspective and axes.