# Inside Koval's 3D Grapher

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

And another:

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\)

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)\).

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 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?

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.

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.

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

## 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:

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.

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:

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:

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

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