Non-Linear Interpolation

Intro
Interpolation is the process of calculating the positions of points at regular intervals between two points, one at a time. For example, a line drawing algorithm takes 2 points as parameters, then it must calculate the exact position of each pixel on the line segment. Such interpolation is called linear interpolation, because the graph of data calculated this way is a straight line. Non-linear interpolation can be used to make a smooth curve between two or more points, though it is slightly more complicated. However, assuming you're decent at math, it should not be a hard concept to understand. I will first discuss linear interpolation, just to make sure that you're familiar with it. Please note that, in this tutorial, I will leave much room for optimization; that is, I will explain the concepts to you, and let you make it faster if you wish.

Linear Interpolation
Suppose you have 2 points, A and B:

To draw a line segment between these 2 points, we simply find whether the line is more horizontal or more vertical, then we calculate how many pixels to move along each axis through each iteration. A simple, slow way to do this in QB is:

SUB LINEX (X1%, Y1%, X2%, Y2%, C%)

DeltaX% = X2% - X1%
DeltaY% = Y2% - Y1%
IF ABS(DeltaX%) > ABS(DeltaY%) THEN
	IncX! = SGN(DeltaX%)
	IncY! = DeltaY% / DeltaX%
	L% = ABS(DeltaX%)
ELSE

	IncX! = DeltaX% / DeltaY%
	IncY! = SGN(DeltaY%)
	L% = ABS(DeltaY%)
END

X! = X1%
Y! = Y1%
FOR A% = 0 TO L%
	PSET (X!, Y!), C%
	X! = X! + IncX!
	Y! = Y! + IncY!
NEXT A%

END SUB

That is obviously quite simple, but illustrates very simple interpolation. Note, also, that the graph formed by this process can be described by a mathematical formula. This is very important, as it allows us to calculate the position of any given point along the line segment. Of course, this formula is:

y = mx + b

And we can change the m and b to be in terms of X1, Y1, X2, and Y2. However, that is irrelevant to this tutorial, and it's rather easy anyway.

Non-Linear Interpolation
As you probably already know, an interpolation function must be mathematically describable, so we must find functions that represent a smooth curve. Three such functions that I know of are cubic, cosine, and parabola. I will only discuss cosine and parabolic interpolation in this tutorial. In general, we'll have to take only part of the function we're using and fit the interpolation formula to that one part.

Cosine Interpolation
As I'm sure you probably know if you're reading this tutorial, the function f(x) = cos x looks like:

However, we want to map only the first half of the cosine function to the space between points A and B:

Of course, our data, A and B, will seldom fit perfectly to the unmodified cosine function, so we take into account the distance between A and B along the X and Y axes. Since we want to calculate the Y coordinate of a point between A and B, we will specify a value, X, between 0 and 1, which will tell how close to A or B the point is. Notice in the picture above where I mapped the cosine function to the points A & B that the cosine function is upside down, assuming B is greater than A (I'll derive the formula with this assumption). Also, the cosine is a periodic function with amplitude 1; that is, the values of the cosine function range from 1 to -1, giving a maximum distance from the X axis of 1 (the amplitude). Of course, this means that the values have a range of 2; however, in order to map the height to the difference in Y values of A & B, we must change it to a range of 1, and we need the values to always be positive. To to this, we can subtract the cosine from 1 and divide it by 2. This will map the height correctly to the data. I will call the distance between A and B along the X axis dx (delta X). Since we're only mapping the first half of the cosine graph, the angle for which to calculate the cosine will be between 0 and pi (180°). Therefore, by multiplying X by pi, we get the right value to get the cosine for. In the following equations, I will assume that A and B are the Y coordinates of the points A and B and that dx is known. In case you want to see exactly what I've been rambling on about:

Now we can develop the formula:

Y = .5 * (1 - cos(X*pi)) * (B - A)

Of course, X can be rewritten as some value, xx/dx, where xx is
between 0 and dx, but that requires a division (ick!). Now we can multiply the 2
quantities above, then factor out A and B:

Y = .5 * (B - A - B * cos(X*pi) + A * cos(X*pi))
Y = .5 * (A * (cos(X*pi) - 1) + B * (1 - cos(X*pi)))

COS() is not the fastest of functions, so, if you wish, you may rewrite the above equation to
calculate the cosine only once:

cosxpi = cos(X * pi)
Y = .5 * (A * (cosxpi - 1) + B * (1 - cosxpi)))

This still isn't an extremely fast routine, though, due to the use of floating point math and three multiplications. It can be optimized by using a look-up table for cosine and changing the multiplication by .5 to SHR. ; )

Parabolic Interpolation
I figured out this method of interpolation by myself, and it's only slightly more complicated than cosine interpolation. I got the idea from my particle engine, because I noticed that, when simulating gravity, the particles followed a parabolic curve, even though I had not specifically "told" it to. The reason it followed a parabolic curve is that I added a value to the Y coordinate of each particle, then increased the amount to add each frame. An over-simplified loop, where y is the particle's Y coordinate, yi is the number to add to y each frame, and i is a "gravity constant" that is added to yi each frame, would look something like:

DO	'main loop
	PSET (x, y), colour	'don't worry about x or colour right now, since a downward force
				'would not affect either of them
	y = y + yi
	yi = yi + i
LOOP

As you probably guessed from the title of this method of interpolation, we will map a parabola to the space between the points A and B. In case you don't remember, a parabola (y = x²) looks something like:

We will actually map two halves of a parabola to the points, but one half will be upside-down. The next part may seem confusing at first, but it's really not bad. What we do now is make an equation for the height, dy, of a curve with width dx. The curve looks like:

Note that this curve is formed using the same method as in my particle engine. Therefore, we can assume that the following is true:

In this equation, a is a constant value that tells us how much to vertically stretch or shrink the curve. Since it's constant, a can be factored out of the summation. So, rearranging the equation to solve for a, which is the unknown variable:

I know this looks nasty, but it can be simplified. This simplification can best be illustrated by a simple diagram of the sum of n as n goes from 1 to, say, 5. Notice that what we're interested in this diagram is the area of the resulting triangle.

Note that I've dotted in the corners of the other spaces which would form a square. Also, the area of the sum of n as n goes from 1 to dx is equal to half the area of the square formed (.5*dx²) plus one half of the length of the square (because the area also includes an extra half of the diagonal). Therefore, we can say:

Now we substitute this into the equation for a:

Here comes the easy part. Since we want to map a parabola to the data, we will use the formula for a parabola:

y = a(x - h)² + k

We don't want to translate the function, so h and k are both zero, making the formula:

y = ax²

Now we just substitute our equation for a into that formula:

But remember that we will actually map 2 half parabolas, so this equation will work only for the first half. Modifying the equation to work for the second half isn't hard, as shown below:

The way it is now, you must check if x less than dx It is possible to combine the two formulas into one, but I will leave that to you (sorry, I have a deadline to make). The interpolation obtained from this method is actually a little more rounded than cosine interpolation, making it a little better looking, in my opinion. Of course, you pay for it with the division and comparison, so this may not be the best method to use. I just found it quite interesting and decided to write about it. The resulting curve will look something like:

Following is a simple routine to use parabolic interpolation between 2 points:

SUB INTERPOLATE (X1, Y1, X2, Y2, C)

DX = X2 - X1
DY = Y2 - Y1
XF = DX * (DX + 1) / 2
IF XF = 0 THEN EXIT SUB
FOR X = 0 TO DX / 2 STEP SGN(DX)
	Y = X * X * DY / XF
	PSET (X1 + X, Y1 + Y), C
	PSET (X1 + DX - X, Y2 - Y), C
NEXT X

END SUB

Conclusion
There are advantages and disadvantages to both methods of non-linear interpolation described here. Which one to use depends both on its application and the programmer's preference. Cosine interpolation has the advantage of having no divisions, yet both methods use floating point math. Parabolic interpolation can be faster for drawing a curve, since you can easily plot 2 points at a time. Of course, you could probably do something similar with cosine interpolation if you wanted. The idea of parabolic interpolation isn't the greatest of ideas, but I think that it's interesting to understand how it works. Of course, only a geek like me could appreciate such a thing. ; ) If you have any questions, comments, corrections, etc., feel free to email me at zippy_@hotmail.com. I hope this all makes sense to you.