```
Jiho Kim
```

Mathematical Graphics

The MSRI Summer School, Reed College

June 18 - July 2, 2005

*I saw an angel in the marble and carved until I set him
free.*

~ Michelangelo ~

Given a function `f(x,y)`

, draw the level set defined by
`f(x,y)=0`

.

We restrict our problem to the following conditions:

`f`

is smooth. (Has second partial derivatives.)- We are drawing on rectangular region of the xy plane, with its sides parallel to an axis. This is not mathematically necessary, but it makes the implementation of the algorithm easier.

The following are examples of functions we might deal with.

`f(x,y)=0` | |

`f(x,y)=a x + b y + c` | |

`f(x,y)=x^2+y^2-r^2` | |

`f(x,y)=y^2-(x-a)(x-b)(x-c)` | |

`f(x,y)=(x^2+y^2)^2-c^2(x^2-y^2)` | |

`f(x,y)=sin(x)+sin(y)` | |

`f(x,y)=sin(x^2+y^2)+sin(x)` |

**Method 1**

Color each pixel one color or another depending on if the functional value is positive or negative. Then, apply edge detection. The following pictures illustrate the result. The green portion of the first image is positive while the black part is negative. Running the image through a edge-detection filter in an open source image manipulation program the GIMP, produced the second image.

**Method 2**

For each pixel, compute the functional value. If the magnitude of that value is less than a given threshhold, paint the pixel.

For `f(x,y)=(x^2+y^2)^2-4(x^2-y^2)`

and a threshold of
1, we get the first picture. In the first picture, we notice that the
'curve' thins out more as it gets farther away from the origin. The
second picture was produced by the same algorithm where the threshold
at given pixel was given by the absolute value of the x coordinate.
This evens out the curve's thickness, but such considerations cannot
be made for an arbitrary function. We also notice as the x-value goes
to 0, there is a gap in the curve since the functional value is
greater than the absolute value of the x-coordinate.

**Can we do better?**

Both methods requires that the functional value of each pixel to be computed. Can we do better than that?

Restrict our functions to be linear. There are many ways for a plane to intersect a rectangular region of the xy-plane.

By computing the functional values at each corner of the rectangular region (seen in blue above), we can figure out which case we must consider. In each case, there is an obvious way to draw a line on the blue square to match the level set of the red surface.

Our algorithm divides the rectangular region into sufficiently small pieces the surface each subregion is closely approximated by a plane.

On the left, we have the true circle. While on the right, we have the approximation gotten by dividing the square along the gray lines (into a 6 by 6 grid). The fact that we have a fixed grid to start with makes this a non-adaptive algorithm.

Of course there is a difference in the two pictures since the surface is far from planar in the subregions. We can compare them better by superimposing these images atop each. We get the following picture.

We notice that the approximation draws a curve inside the true curve. The next section discusses why that is.

The applet first shows the lemniscate given by `(x^2+y^2)^2-3(x^2-y^2)=0`

. Clicking on the graph will flip through approximations of the lemniscate with a given grid size. With a 5x5 grid, we can see that the algorithm completely fails to pick up the saddle at the origin. Even with a 32x32 grid, we can see artifacts of the grid near the origin. The adaptive algorithm divides the regions near the origin more to remove the zig-zig in the curve that we observe with the 32x32 grid.

0. Initialize the stack with the whole region to be considered.

1. `While`

the region stack is not empty, pop the top one off
the stack and repeat steps 2-3.

2. `If`

the region is *small enough* then apply the
linear case and draw the segment.

3. `Else`

, if the region is bigger than some threshold,
subdivide the region into smaller regions,
push the subregions onto the stack. If the region is smaller than
the threshold, discard it.

The following applet how a rectangular region can be subdivided to focus on certain points. This is not exactly what happens with the algorithm for drawing contour curves, but idea is the same.

But we still need to talk about what is small enough.

Let us assume that the function can be approximated by a second degree Taylor polynomial. This requires that the function has second derivatives. We know that if the second degree coefficients of the polynomial are negligible (all close to 0) then the approximation is almost a linear function, i.e. a plane. So we must reduce the size of the region until the sum of the absolute value of degree-2 coefficients in the Taylor polynomial is below some threshold. This number equals the norm of the Hessian matrix.

Assuming that the surface defined on a small region is almost planar, then we must consider the magnitude of the gradient. The magnitude of the gradient is inverse proportional to the error introduced by some fixed vertical perturbation of the function.

You can see in the applet above that closer the slope is to 0, the greater the horizontal error there is given a fixed vertical perturbation (of 1 unit, in this case).

So must also take into account the inverse of the magnitude of the gradient.

In practice, we consider the product of these factors--the norm of the Hessian matrix and the magnitude of the gradient. In this way, we can get a nonnegative coefficient for each point in a region. Given a region, if the supremum of this product over the region is less than some threshold, we consider the region to be "small enough".

`Lemniscate`

`Elliptic Curve`

`Grid`

`Silly`