Lagrange interpolation is a way of crafting a function from a set of data points..

In the past I’ve seen reference to Lagrange interpolation in relation to audio programming like, for helping make a soft knee for a limiter, but it can be used wherever you need to make a function from some data points.

What’s It Do?

Lagrange interpolation is a way of crafting a function from a set of data pairs. The resulting function passes through all the data points you give it (like a Catmull-Rom spline does), so can be used to find a function to interpolate between data sets.

You can’t give two value pairs that have the same x value, but the data points don’t have to be evenly spaced.

Also, if you give data points, you’ll get out a function that is a degree polynomial. So, if you interpolate two data points, you’ll get a degree 1 polynomial (a line). If you interpolate three data points, you’ll get a degree 2 polynomial (a quadratic).

The function will be quite messy, but you can use algebra, or wolframalpha.com (or the like) to simplify it for you to a simpler equation.

Lagrange interpolation is subject to Runge’s Phenomenon, so the more data points you have, the more the interpolation tends to get “squirly” near the edges and shoot off up high or down low, instead of smoothly interpolating between data values.

How’s It Do It?

Well, to make any kind of curve from data points, if we want the curve to pass through those data points, one way would be to come up with a set of functions to multiply each data point by.

Each function must evaluate to 1 when the curve is at that control point, it should be zero when the curve is at any other control point. Between control points, the function can take any value, but if you make it continuous / smooth, the curve will be continuous and smooth, so that’s usually what is desired.

When we have those functions, to get a point on the curve we just multiply each control point by it’s corresponding function (called a basis function), and we sum up the results.

The pseudocode below is how this works and is the basic functionality of most common curve types:

// The basic way to evaluate most any type of curve
float PointOnCurve (float t, float *controlPoints, int numControlPoints)
{
    float value = 0.0f;

    for (int i = 0; i < numControlPoints; ++i)
        value += controlPoints[i] * ControlPointFunction(i, t);

    return value;
}

float ControlPointFunction (int i, float t)
{
  // return the ith control point function evaluated at time t.
  // aka return f(t) for the ith control point.
}

What makes Lagrange interpolation different than other curve types is the basis functions it uses.

The Math

If you aren’t used to seeing a capital pi, or a laplacian style cursive l in equations, it’s about to get a bit mathy!

If you feel like skipping to the next section, I don’t blame you, but if you are feeling brave, you should try and follow along, because I’m going to slowly walk through each and every symbol to help explain what’s going on and why.

Let’s say that you are want to be able to interpolate between data points:

The formula for calculating a Lagrange interpolated value is this:

The capital sigma () just means that we are going to loop a variable j from 0 to k (including k), and we are going to sum up the total of everything on the right for all values of j. When you see a capital sigma, think sum (note they both start with an s).

The next thing after the sigma is . That is just the y value from our jth control point. That is essentially controlPoints[j].y.

After that comes the last part . That is just the function for the jth control point that we multiply the control point by (aka the basis function), evaluated for the specific value x.

Since there is no operator between this function and the control point, that means we multiply them together. So yeah… that crazy math just says “multiply each control point by it’s basis function, and sum up the results”, just like our pseudo code above does!

The second equation we need to look at is the definition of the basis functions for each control point. Here is the formula that describes the jth basis function, for the jth control point:

First is the capital pi . This means that we are going to do a loop, but instead of adding the results of the loop, we are going to multiply them together. Where a capital sigma means sum, capital pi means product.

The notation for product is a bit different here than in the sigma though which may be a bit tricky to read at first. Instead of explicitly saying that should go from 0 to , the notation $latex 0\le m\le k\\$ says that implicitly. That same notation can be used with sigma, or the more explicit style notation could be used with pi.

The pi also has this notation next to it . That means that the case where equals should be skipped.

Finally, on to the second part: . This part is pretty easy to read. is the parameter to the function of course, is just controlPoints[m].x where is the index variable of our product loop (), and is just controlPoints[j].x where is the index variable of our summation loop ().

Let’s say that was 2 because we had 3 data pairs. Our three basis functions would be:

Which means that our final Lagrange interpolation function would be:

That is quite a mouth full, but hopefully you understand how we came up with that!

is just controlPoints[i].x and is just controlPoints[i].y.

Math Intuition

The intuition here is that we need to come up with a set of functions to multiply each control point by, such that when the function’s x value is at the control point’s x value, the function should evaluate to 1. When the function’s x value is at a different control points x value, the function should evaluate to 0. The rest of the time, the function can evaluate to whatever it wants, although again, having it have smooth values is nice to making a good curve.

So the first problem is, how do we make a function evaluate to 0 when x is at a different control point?

The easy way would be to multiply a bunch of terms together of this form , but make sure and not include the x of the actual control point that we are multiplying against.

That is exactly what it does with the numerator in the product notation of the basis function.

Note that is the index of the current control point that we are calculating the basis function for. All values of x, that isn’t the x value of a control point will evaluate to non zero.

The denominator value is there so that when x is the value of the control point that we care about, that the function will evaluate to 1.

It does this by figuring out what the value of the numerator will be when x is at the control point, and then makes that be the value that it divides by, so that it’s 1 at that x value.

Not too much to it. Pretty simple stuff, but powerful as well!

Extending to 2D and Beyond

Lagrange interpolation is a one dimensional interpolation scheme, meaning that if you have data points of the form (x,y), it can give you an interpolated y value based on an x value you give it. The interpolation it does can never give two different y values for the same x.

If you want to extend this technique to interpolating a curve through two dimensional data points, or even higher, you need to do interpolation independently for each axis and use a “parametric” value for that axis.

For instance, if you needed to interpolate a curve through 3 dimensional points, you would have data points like this:

X Points =
Y Points =
Z Points =

And then you would interpolate on each axis by the t value to get your X, Y and Z axis values. This should look familiar, because this is how higher dimensional Bezier curves work; you evaluate them per axis based on a parametric value per axis (s,t,u,etc).

You could use the same t values on each axis, or they could be completely independent. You don’t even need to have the same number of points for each axis!

You might wonder how this differs from the standard interpolation in the 2D case. Check the demos in the link section below to really get a grasp of the difference, but in essence, with standard (1D) interpolation, you can never have two x values that evaluate to 2 different y values. Extending it like the above into two dimensions by parameterizing each axis lets you get around that limitation and you can make true 2d shapes.

Lastly, it is possible to make Lagrange interpolated surfaces! I won’t go into the details (perhaps a future post!), but if you know how to make a bezier rectangle by doing a tensor product (basically having X axis Bezier curves, multiplied by Y axis Bezier curves), you can accomplish a Lagrange surface in a really similar way.

Sample Code

This sample code is written for readability, but could easily be optimized for faster execution. Also, from what I hear, the second form of Barycentric Lagrange Interpolation is touted as the fastest form of Lagrange interpolation, since many values can be pre-calculated and re-used for different values of x.

#include <stdio.h>
#include <vector>

struct SPoint
{
    float x;
    float y;
};

typedef std::vector<SPoint> TPointList;

void WaitForEnter ()
{
    printf("Press Enter to quit");
    fflush(stdin);
    getchar();
}

// calculates the lagrange basis function y value for control point "controlPointIndex" at x value "x"
float LagrangeBasis (const TPointList& pointList, size_t controlPointIndex, float x)
{
    // this is the pi "inner loop" multiplication work
    float value = 1.0f;
    for (size_t i = 0, c = pointList.size(); i < c; ++i) {
        if (i != controlPointIndex)
            value *= (x - pointList[i].x) / (pointList[controlPointIndex].x - pointList[i].x);
    }
    return value;
}

// returns a value at x, using lagrange interpolation over the specified list of (x,y) pairs
float LagrangeInterpolate (const TPointList& pointList, float x)
{
    // this is the sigma "outer loop" summation work
    float sum = 0.0f;
    for (size_t controlPointIndex = 0, c = pointList.size(); controlPointIndex < c; ++controlPointIndex)
        sum += pointList[controlPointIndex].y * LagrangeBasis(pointList, controlPointIndex, x);
    return sum;
}

int main (int argc, char **argv)
{
    // show some 1d interpolated values
    // note that the points don't need to be sorted on X, but it makes for easier to read examples
    {
        // (x,y) pairs
        const TPointList points =
        {
            { 0.0f, 1.1f },
            { 1.6f, 8.3f },
            { 2.3f, 6.5f },
            { 3.5f, 4.7f },
            { 4.3f, 3.1f },
            { 5.9f, 7.5f },
            { 6.8f, 0.0f }
        };

        // show values interpolated from x = 0, to x = max x
        printf("1d interpolated values.  y = L(t)n");
        const float c_numPoints = 10;
        for (int i = 0; i < c_numPoints; ++i)
        {
            float percent = ((float)i) / (float(c_numPoints - 1));
            float x = points.back().x * percent;
            float y = LagrangeInterpolate(points, x);
            printf("  (%0.2f, %0.2f)n", x, y);
        }
        printf("n");
    }

    // show some 2d interpolated values
    // also note that x and y don't have to have matching t values!
    {
        // (t, x) pairs
        const TPointList pointsX =
        {
            { 0.0, 0.0f},
            { 1.0, 1.6f},
            { 2.0, 2.3f},
            { 3.0, 3.5f},
            { 4.0, 4.3f},
            { 5.0, 5.9f},
            { 6.0, 6.8f}

        };
        // (t, y) pairs
        const TPointList pointsY =
        {
            { 0.0f, 1.1f },
            { 1.0f, 8.3f },
            { 2.0f, 6.5f },
            { 3.0f, 4.7f },
            { 4.0f, 3.1f },
            { 5.0f, 7.5f },
            { 6.0f, 0.0f }
        };

        // show values interpolated from t = 0, to t = max t, on each axis
        printf("2d interpolated values.  x = L(t_x), y = L(t_y)n");
        const float c_numPoints = 10;
        for (int i = 0; i < c_numPoints; ++i)
        {
            float percent = ((float)i) / (float(c_numPoints - 1));

            // calculate x
            float tx = pointsX.back().x * percent;
            float x = LagrangeInterpolate(pointsX, tx);

            // calculate y
            float ty = pointsY.back().x * percent;
            float y = LagrangeInterpolate(pointsY, ty);

            printf("  (%0.2f, %0.2f)n", x, y);
        }
        printf("n");
    }

    WaitForEnter();
    return 0;
}

And here’s the programs output:

Final Notes

Now that you know how to do all this stuff I wanted to share a couple more pieces of info.

Firstly, it’s kind of weird to call this “Lagrange Interpolation”. A better term is to call this the “Lagrange Form of Polynomial Interpolation”. The reason for that is that if you have some number of data points, there exists only one unique minimal order polynomial (lowest degree of x possible) that fits those points. That is due to the “unisolvence theorem” that you can read more about here: Wikipedia: Polynomial interpolation.

What that means is that if you were to use a different type of polynomial interpolation – such as newton interpolation – the result you get out is algebraically equivalent to the one you’d get from this Lagrange form. There are pros and cons to using different forms of polynomials, but that’s out of the scope of this post so go read about them if you are interested!

Speaking of that, even though this sample code is focused on interpolation using the Lagrange form, this technique is really great at being able to just come up with some simpler f(x) function that passes through specific data points. In this way, you can kind of “bake out” a custom f(x) function to do interpolation for specific values, that doesn’t need all the moving parts of the Lagrange form. For example, if you make the formula for lagrange interpolation of 3 specific value pairs and then simplify, will get out a simple quadratic function in the form of !

Links

Here are some interactive demos I made to let you play with Lagrange interpolation to get a feel for how it works, and it’s strengths and weaknesses:
One Dimensional Lagrange Interpolation
Two Dimensional Lagrange Interpolation

I also found these links really helpful in finally understanding this topic:
Lagrange Interpolation
Lagrange’s Interpolation Formula

Want to follow the rabbit hole a little deeper? Check out how sinc interpolation relates to the Lagrange form!
The ryg blog: sinc and Polynomial interpolation

The usual equation for blending two pixels with alpha is to use a source factor of “source alpha” and a destination factor of “one minus source alpha”.

That results in this equation:

In other words, it’s just a linear interpolation between and , using the alpha of the pixel you are writing to determine the weighting for the lerp.

Just to be clear, is the pixel in your frame buffer (ie the final result that shows up on your monitor) and is the new pixel you are trying to write or combine with the output image.

If you use pre-multiplied alpha, all that means is that is already multipled by which results in slightly less math:

Less math for the same results is always a good thing due to increased efficiency, but this also results in a higher quality results in the case of using mips. For more info on that check out this link: NVIDIA GameWorks – Alpha Blending: To Pre or Not To Pre

To achieve this in OpenGL or DirectX with pre-multiplied alpha textures, you just use a source factor of “one” and leave the destination factor at “one minus source alpha”.

If you are using alpha to write to a render target that also has an alpha channel, the math improvement is even better. Here is the equation with regular (post-multiplied) alpha in that scenario for combining the RGB portion of the pixels:

When working with pre-multiplied alpha it just goes back to the lerp:

When blending to a target with alpha, the equation to combine alpha is the old familiar lerp:

That makes the lerp form of color combining even nicer since it means the same math for all color channels of the pixel:

Confused about why alpha combining using lerp makes sense? It makes most sense to me when thinking about it like this… if you had something that was half opaque(Out.A = 0.5), then you looked at that through something that was 3/4 opaque (In.A = 0.75), only 25% of the light would get past the top layer (0.25), to reach the bottom layer. Only 50% of the light that reached the bottom layer gets through, so we cut that 25% in half to get 12.5% (0.125). Since alpha really means “opaqueness” (1.0 means no transparency), that means that the combined alpha is 1 – 0.125, or 0.875.

If you plug the numbers into the above equation you get the same result:

This compute graphics stack exchange question has some more great info on what premultiplied alpha can give you. It’s much more than I wrote here:
Does premultiplied alpha give order independent transparency?

Here’s a fun question… does alpha represent transparency, or does it represent how much of a pixel is covered by an opaque object? To find out the answer give this a read!
Jounral of Computer Graphics Techniques: Interpreting Alpha

In a previous post, I showed Four Ways to Calculate Sine Without Trig. While reading up on rational bezier curves, I came across a 5th way!

A rational bezier curve post will be coming in the (hopefully near) future, so I won’t go into too many details of that, but the important thing to know is that with a rational 2d quadratic bezier curve, you can EXACTLY represent any conic section, including circle arcs. With a regular (integral) bezier curve, you can’t.

You can in fact see that in action in an interactive HTML5 demo i made here: Rational Quadratic Bezier Curve. Note that on my main page, i have links to other 1d/2d bezier/trig rational/integral interactive curve demos. go check them out if you are interested in seeing more: http://demofox.og.

After I made the 2d quadratic bezier demo I realized something… The x and y positions of that curve are calculated independently of each other, and if you wanted to draw a circle in a program, you’d calculate the x and y positions independently, using sine for the y axis value and cosine for the x axis value. That means that if rational curves can – as people say – perfectly represent conic sections, that the two 1d curves that control each axis of that circle arc curve must be exact representations of sine and cosine, at least for the 90 degrees of that arc.

It turns out that this is true enough. I am seeing some variation, but have an open question on math stack exchange to try to get to the bottom of that.

Let’s do some math to come up with our curve based sine equation!

Math Derivation – Skip If You Want To

A rational Bezier equation is defined like this, which is basically just one Bezier curve divided by another:

So, here is the equation for a rational quadratic bezier curve (n = 2):

are the control points and are the weightings associated with those control points.

For the first 90 degrees of sine, .

arcAngle is 90 degrees in our case, so or or .

If we plug those values in and simplify, and treat it as an explicit (1d) equation of , instead of , we come up with the equation in the next section.

Final Equation

or, so you can copy paste it:
y=(0.70710678118*2x(1-x) + x^2) / ((1-x)^2 + 0.70710678118*2x(1-x) + x^2)

That gives us the first quadrant (first 90 degrees) of sine. To get the second quadrant, we just horizontally flip the first quadrant. To get the third quadrant, we vertically flip the first quadrant. To get the fourth quadrant, we vertically and horizontally flip the first quadrant.

Equation In Action

Here’s a screenshot from a shadertoy where I implemented this. red = true sine value, green = the value made with the curve, yellow = where they overlap and are equal. Any place you see red or green peaking out means they aren’t quite equal. I sort of expected it to be exactly dead on correct, but like I said, I have an open math stack exchange question to find out why it isn’t and will update this with my findings if they are relevant or interesting. It still is pretty darn good though.

Also, this formula may not be the most efficient way to calculate sine without trig that I’ve shown, but it is interesting, and that’s why I wanted to show it (:

You can see this in action at Shadertoy: Sin Without Trig IV

Take a look at the picture below. Can you calculate the values at points A,B and C?

There’s a technique for doing this called bilinear interpolation which extends the idea of linear interpolation into two dimensions.

To calculate bilinear interpolation, you just do linear interpolation on one axis, and then the other.

Check out these examples:

Point A
Point A has a coordinate of (0.2,0.8). Let’s start with the X axis.

Firstly, we’ll do interpolation across the top of the square on the X axis. That means we need to go 0.2 (20%) of the way from 7 to 0. That 0.2 is our x axis coordinate. To calculate that we do this:
7 + (0-7) * 0.2 = 5.6

Next, we’ll do interpolation across the bottom of the square on the X axis. So similar to above, we’ll go 0.2 (20%) of the way from 3 to 5. To calculate that we do this:
3 + (5-3) * 0.2 = 3.4

Now that we have interpolated across our X axis we need to interpolate across our Y axis. We are now going to go 0.8 (80%) of the way from our bottom value (the value of Y = 0), to our top value (the value of Y = 1). The bottom value is the bottom interpolation we did (3.4) and the top value is the top interpolation we did (5.6).

So, we need to go 0.8 (80%) of the way from 3.4 to 5.6. To calculate that we do this:
3.4 + (5.6 – 3.4) * 0.8 = 5.16

So, the value at point A is 5.16. The image below should help explain visually the process we went through.

Point B

Point B has a coordinates of (1.0,0.5).

Let’s start with the X axis interpolation again (although you could easily start with Y first if you wanted to! You’d end up with the same result).

Doing the X axis interpolation across the top, we need to go 1.0 (100%) of the way from 7 to 0. I’m sure you can guess what the answer is but here it is calculated:
7 + (0-7) * 1.0 = 0

Similarly, let’s do the X axis interpolation across the bottom. We need to go 1.0 (100%) of the way from 3 to 5.
3 + (5-3) * 1.0 = 5

Next up comes the Y axis interpolation. We need to go 0.5 (50%) of the way from 5 to 0.
5 + (0-5) * 0.5 = 2.5

There we go, the value at point B is 2.5. That should make sense too by looking at the diagram. It’s basically 1d linear interpolation between 5 and 0, and is half way between them, so intuitively, the answer 2.5 should make sense.

Point C

Point C has a coordinate of (0.8,0.2).

Once again, let’s do the X axis interpolation across the top of the box. We need to go 0.8 (80%) of the way from 7 to 0.
7 + (0-7) * 0.8 = 1.4

Then, we need to do the x axis interpolation across the bottom of the box. We need to go 0.8 (80%) of the way from 3 to 5.
3 + (5-3) * 0.8 = 4.6

Lastly, the y axis interpolation. We need to go 0.2 (20%) from 4.6 to 1.4.
4.6 + (1.4-4.6) * 0.2 = 3.96

The value of point C is 3.96

Bilinear Filtering

While bilinear interpolation is useful if you have a grid of data that you want to interpolate values within (such as a height field that represents terrain?!), where this really shines in game development is in bilinear texture filtering.

Basically, when you have a texture coordinate that doesn’t perfectly line up with the center of a pixel in the texture (because there’s a remainder), it uses the fractional part of the coordinate within that pixel to do bilinear interpolation between the 4 pixels involved to come up with the final pixel color. It does bilinear interpolation on each of the channels: Red, Green, Blue and Alpha.

The end result is pretty good! Check out the image below to see it in action on a texture that we are zoomed in too far on. The left side uses bilinear filtering, while the right side does not, and just shows you the nearest pixel.

Consequently, bilinear texture filtering has to do 4 pixel reads to be able to interpolate between them, instead of a single pixel read like when doing nearest pixel sampling. More pixel reads = more overhead, and not as cheap (computationally) to use.

Links

For some folks who are reading this, this info is pretty basic and you might be wondering why i bothered writing about it. Look to the next post for something kind of bizarre regarding bilinear filtering. This was needed as a stepping stone to the next thing I want to show you guys 😛

Wikipedia: Linear Interpolation

Wikipedia: Bilinear Interpolation

Wikipedia: Bilinear Filtering

The Dark Side of Bilinear Filtering

Check out these links for some deeper info about some shortcomings (and some workarounds) of hardware based bilinear sampling:
iq: hardware interpolation
iq: improved texture interpolation