An Adaptive Algorithm for Drawing Level Curves

Jiho Kim
Mathematical Graphics
The MSRI Summer School, Reed College
June 18 - July 2, 2005

GCalc: Java Graphing Calculator

Motivation

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

~ Michelangelo ~


The Problem

Given a function f(x,y), draw the level set defined by f(x,y)=0.

We restrict our problem to the following conditions:


Examples

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)

Naive Approaches

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.

Before edge detection After edge detection

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.

Before edge detection After edge detection

Can we do better?

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

Simplified Problem

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.

Non-adaptive algorithm

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

true circle approximation with 6x6 grid

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.

composite of true circle with approximation

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

What goes wrong...

More things that go wrong...

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.

Adaptive Algorithm

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.

When is a region 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".

More Pictures

Lemniscate

More Pictures

Elliptic Curve

More Pictures

Grid


More Pictures

Silly