Finite differences are numerical methods for approximating function derivatives – otherwise known as the slope of a function at a specific point on the graph. This can be helpful if it’s undesirable or impossible to calculate the actual derivative of a specific function.

This post talks about three methods: central difference, backwards difference and forward difference. They are all based on evaluating the function at two points and using the slope between those points as the derivative estimate.

The distance between those sample points is called an epsilon, and the smaller it is, the more accurate the approximation is in theory. In practice, extremely small values (like FLT_MIN) may hit numerical problems due to floating points number usage, and also you could hit performance problems due to using floating point denormals. Check the links section at the bottom for more info on denormals.

Central Difference

The central difference is the most accurate technique of the three. You can find information about comparitive accuracy of these three techniques in the links section at the end. In practical terms, this may also be the slowest method too – or the most computationally expensive – which i’ll explain further down.

If you want to know the derivative of some function at a specific value of x, you pick an epsilon e and then you calculate both and . You subtract the first one from the second and divide by 2*e to get an approximated slope of the function at the specific value of x.

Remembering that the slope is just rise over run, and that the derivative at a point on a function is just the slope of the function at that point, this should hopefully make sense and be pretty intuitive why it works.

The resulting equation looks like:

This process is visualized below. The black line is the actual slope at 0.4 and the orange line is the estimated slope. The orange dots are the sample points taken. The epsilon in this case is 0.2.

Interestingly, when dealing with quadratic (or linear) functions, the central difference method will give you the correct result. The picture above uses a quadratic function, so you can see no matter what value of e we use, it will always be parallel to the actual slope at that point. For cubic and higher functions, that won’t always be true.

Backward Difference

The backward difference works just like the central difference except uses different sample points. It evaluates and , subtracts the 1st one from the second one and divides the result by e.

The resulting equation looks like this:

A neat property shared by both this and the forward difference is that many times you are already going to be evaluating f(x) for other uses, so in practice this will just mean that you only have to evaluate f(x-e), and will already have the f(x) value. That can make it more efficient than the central difference method, but it can be less precise.

Also, if you are walking down a function (say, rendering a Bezier curve, and wanting the slope at each point to do something with), you may very well be able to use the f(x) of the previous point as your f(x-e) function, which means that you could possibly calculate the backwards difference by using the previous point, instead of evaluating the function extra times in your loop!

Check out the image below to see how different values of e result in different quality approximations. The smaller the epsilon value, the more accurate the result. An infinitely small epsilon would give the exact right answer.

Forward Difference

The forward difference is just like the backwards difference but it evaluates forward instead of backwards.

The equation looks like this:

Below you can see it visually. Note again that smaller values of e make the estimation closer to correct.

On the GPU

If you’ve ever encountered the glsl functions dFdx and dFdy and wondered how they work, they actually use these same techniques.

Shaders run in groups, and using dFdx, the shader just looks to it’s neighbor for the value that was passed to it’s dFdx, then using “local differencing” (per the docs), gives each shader the derivative it was able to calculate.

Links

GLSL: dFdx, dFdy
Wikipediate: Finite Difference – The wikipedia page talks about more details, including how to calculate 2nd derivatives and higher!
Floating Point Denormals, Insignificant But Controversial
Comparing Methods of First Derivative Approximation Forward, Backward and Central Divided Difference

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