Bicubic Interpolation?

Question

I looked through the internet, and in terms of Bicubic Interpolation, I can't find a simple equation for it. Wikipedia's page on the subject wasn't very helpful, so is there any easy method to learning how Bicubic Interpolation works and how to implement it? I'm using it to generate Perlin Noise, but using bilinear interpolation is way to choppy for my needs (I already tried it).

If anyone can point me in the right direction by either a good website or just an answer, I would greatly appreciate it. (I'm using C# by the way)

Answer 1

Using this (Thanks to Ahmet Kakıcı who found this), I figured out how to add Bicubic Interpolation. For those also looking for the answer, here is what I used:

private float CubicPolate( float v0, float v1, float v2, float v3, float fracy ) {
    float A = (v3-v2)-(v0-v1);
    float B = (v0-v1)-A;
    float C = v2-v0;
    float D = v1;

    return A*Mathf.Pow(fracy,3)+B*Mathf.Pow(fracy,2)+C*fracy+D;
}

In order to get 2D Interpolation, I first got the x, then interpolated the y. Eg.

float x1 = CubicPolate( ndata[0,0], ndata[1,0], ndata[2,0], ndata[3,0], fracx );
float x2 = CubicPolate( ndata[0,1], ndata[1,1], ndata[2,1], ndata[3,1], fracx );
float x3 = CubicPolate( ndata[0,2], ndata[1,2], ndata[2,2], ndata[3,2], fracx );
float x4 = CubicPolate( ndata[0,3], ndata[1,3], ndata[2,3], ndata[3,3], fracx );

float y1 = CubicPolate( x1, x2, x3, x4, fracy );

Where ndata is defined as:

float[,] ndata = new float[4,4];
for( int X = 0; X < 4; X++ )
    for( int Y = 0; Y < 4; Y++ )
        //Smoothing done by averaging the general area around the coords.
        ndata[X,Y] = SmoothedNoise( intx+(X-1), inty+(Y-1) );

(intx and inty are the floored values of the requested coordinates. fracx and fracy are the fractional parts of the inputted coordinates, to be x-intx, and y-inty, respectively)

Answer 2

Note: The Bicubic formula given by Nicholas above (as the answer) has an error in it. by interpolating a sinusoid, I was able to see it was not correct.

The correct formula is:

float A = 0.5f * (v3 - v0) + 1.5f * (v1 - v2);
float B = 0.5f * (v0 + v2) - v1 - A;
float C = 0.5f * (v2 - v0);
float D = v1;

For the derivation, see https://www.paulinternet.nl/?page=bicubic

Answer 3

I'm a bit confused on the third degree polynomial used.

Yes it gives the correct values in 0 and 1, but the derivatives of neighbouring cells does not fit, as far as I can calculate. If the grid-data is linear, it does not even return a line....

And it is not point symmetric in x=0.5

The polynomial that fits in 0 and 1 AND also have the same derivatives for neighbouring cells, and thus is smooth, is (almost) as easy to calculate.

(and it reduces to linear form if that fits the data)

In the labelling p and m is a short hand for plus and minus e.g. vm1 is v(-1)

    //Bicubic convolution algorithm, cubic Hermite spline
    static double CubicPolateConv
            (double vm1, double v0, double vp1, double vp2, double frac) {
        //The polynomial of degree 3 where P(x)=f(x) for x in {0,1}
        //and P'(1) in one cell matches P'(0) in the next, gives a continuous smooth curve.
        //And we also wants it to reduce nicely to a line, if that matches the data
        //P(x)=Dx³+Cx²+Bx+A=((Dx+C)x+B)x+A
        //P'(x)=3Dx²+2Cx+B
        //P(0)=A       =v0
        //P(1)=D+C+B+A =Vp1
        //P'(0)=B      =(vp1-vm1)/2
        //P'(1)=3D+2C+B=(vp2-v0 )/2
        //Subtracting expressions for A and B from D+C+B+A  
        //D+C =vp1-B-A = (vp1+vm1)/2 - v0
        //Subtracting that twice and a B from the P'(1)
        //D=(vp2-v0)/2 - 2(D+C) -B =(vp2-v0)/2 - (Vp1+vm1-2v0) - (vp1-vm1)/2
        // = 3(v0-vp1)/2 + (vp2-vm1)/2
        //C=(D+C)-D = (vp1+vm1)/2 - v0 - (3(v0-vp1)/2 + (vp2-vm1)/2)
        // = vm1 + 2vp1 - (5v0+vp2)/2;
        //It is quite easy to calculate P(½)
        //P(½)=D/8+C/4+B/2+A = (9*(v0+vp1)-(vm1+vp2))/16
        //i.e. symmetric in its uses, and a mean of closest adjusted by mean of next ones

        double B = (vp1 - vm1) / 2;
        double DpC =(vp1 -v0) -B; //D+C+B+A - A - B
        double D = (vp2 - v0) / 2 - 2 * DpC - B;
        double C = DpC - D;
        //double C = vm1 + 2 * vp1 - (5 * v0 + vp2) / 2;
        //double D = (3*(v0 - vp1) + (vp2 - vm1)) / 2;

        return ((D * frac + C) * frac + B) * frac + A;
    }

Inspired by the comment of ManuTOO below I made it as fourth order too, with an optional parameter, you can set/calculate as you please, without breaking the smoothness of the curve. It is basically the same with an added term all the way in the calculations. And If you set E to Zero, it will be identical to the above. (Obviously if the data is actually on a line your calculation of E must assure that E is zero to get a linear output)

        //The polynomial of degree 4 where P(x)=f(x) for x in {0,1}
        //and P'(1) in one cell matches P'(0) in the next, gives a continuous smooth curve.
        //And we also wants the it to reduce nicely to a line, if that matches the data
        //With high order quotient weight of your choice....
        //P(x)=Ex⁴+Dx³+Cx²+Bx+A=((((Ex+D)x+C)x+B)x+A
        //P'(x)=4Ex³+3Dx²+2Cx+B
        //P(0)=A          =v0
        //P(1)=E+D+C+B+A  =Vp1
        //P'(0)=B         =(vp1-vm1)/2
        //P'(1)=4E+3D+2C+B=(vp2-v0 )/2
        //Subtracting Expressions for A, B and E from the E+D+C+B+A 
        //D+C =vp1-B-A -E = (vp1+vm1)/2 - v0 -E
        //Subtracting that twice, a B and 4E from the P'(1)
        //D=(vp2-v0)/2 - 2(D+C) -B -4E =(vp2-v0)/2 - (Vp1+vm1-2v0-2E) - (vp1-vm1)/2 -4E
        // = 3(v0-vp1)/2 + (vp2-vm1)/2 -2E
        //C=(D+C)-(D) = (vp1+vm1)/2 - v0 -E - (3(v0-vp1)/2 + (vp2-vm1)/2 -2E)
        // = vm1 + 2vp1 - (5v0+vp2)/2 +E;

        double E = ????????; //Your feed.... If Zero, cubic, see below
        double B = (vp1 - vm1) / 2;
        double DpC =(vp1 -v0) -B -E; //E+D+C+B+A - A - B -E
        double D = (vp2 - v0) / 2 - 2 * DpC - B - 4 * E;
        double C = DpC - D;

        return (((E * frac + D) * frac + C) * frac + B) * frac + v0;

ADD: Quoted from suggestion by @ManTOO below:

         double E = (v0 - vm1 + vp1 - vp2) * m_BicubicSharpness;

With m_BicubicSharpness at 1.5, it's very close of Photoshop's Bicubic Sharper ; personally, I set it to 1.75 for a bit of extra kick.

Note that if data is on a line, this suggestion reduces nicely to zero

Answer 4

Took Eske Rahn answer and made a single call (note, the code below uses matrix dimensions convention of (j, i) rather than image of (x, y) but that shouldn't matter for interpolation sake):

/// <summary>
/// Holds extension methods.
/// </summary>
public static class Extension
{
    /// <summary>
    /// Performs a bicubic interpolation over the given matrix to produce a
    /// [<paramref name="outHeight"/>, <paramref name="outWidth"/>] matrix.
    /// </summary>
    /// <param name="data">
    /// The matrix to interpolate over.
    /// </param>
    /// <param name="outWidth">
    /// The width of the output matrix.
    /// </param>
    /// <param name="outHeight">
    /// The height of the output matrix.
    /// </param>
    /// <returns>
    /// The interpolated matrix.
    /// </returns>
    /// <remarks>
    /// Note, dimensions of the input and output matrices are in
    /// conventional matrix order, like [matrix_height, matrix_width],
    /// not typical image order, like [image_width, image_height]. This
    /// shouldn't effect the interpolation but you must be aware of it
    /// if you are working with imagery.
    /// </remarks>
    public static float[,] BicubicInterpolation(
        this float[,] data, 
        int outWidth, 
        int outHeight)
    {
        if (outWidth < 1 || outHeight < 1)
        {
            throw new ArgumentException(
                "BicubicInterpolation: Expected output size to be " +
                $"[1, 1] or greater, got [{outHeight}, {outWidth}].");
        }

        // props to https://stackoverflow.com/a/20924576/240845 for getting me started
        float InterpolateCubic(float v0, float v1, float v2, float v3, float fraction)
        {
            float p = (v3 - v2) - (v0 - v1);
            float q = (v0 - v1) - p;
            float r = v2 - v0;

            return (fraction * ((fraction * ((fraction * p) + q)) + r)) + v1;
        }

        // around 6000 gives fastest results on my computer.
        int rowsPerChunk = 6000 / outWidth; 
        if (rowsPerChunk == 0)
        {
            rowsPerChunk = 1;
        }

        int chunkCount = (outHeight / rowsPerChunk) 
                         + (outHeight % rowsPerChunk != 0 ? 1 : 0);

        var width = data.GetLength(1);
        var height = data.GetLength(0);
        var ret = new float[outHeight, outWidth];

        Parallel.For(0, chunkCount, (chunkNumber) =>
        {
            int jStart = chunkNumber * rowsPerChunk;
            int jStop = jStart + rowsPerChunk;
            if (jStop > outHeight)
            {
                jStop = outHeight;
            }

            for (int j = jStart; j < jStop; ++j)
            {
                float jLocationFraction = j / (float)outHeight;
                var jFloatPosition = height * jLocationFraction;
                var j2 = (int)jFloatPosition;
                var jFraction = jFloatPosition - j2;
                var j1 = j2 > 0 ? j2 - 1 : j2;
                var j3 = j2 < height - 1 ? j2 + 1 : j2;
                var j4 = j3 < height - 1 ? j3 + 1 : j3;
                for (int i = 0; i < outWidth; ++i)
                {
                    float iLocationFraction = i / (float)outWidth;
                    var iFloatPosition = width * iLocationFraction;
                    var i2 = (int)iFloatPosition;
                    var iFraction = iFloatPosition - i2;
                    var i1 = i2 > 0 ? i2 - 1 : i2;
                    var i3 = i2 < width - 1 ? i2 + 1 : i2;
                    var i4 = i3 < width - 1 ? i3 + 1 : i3;
                    float jValue1 = InterpolateCubic(
                        data[j1, i1], data[j1, i2], data[j1, i3], data[j1, i4], iFraction);
                    float jValue2 = InterpolateCubic(
                        data[j2, i1], data[j2, i2], data[j2, i3], data[j2, i4], iFraction);
                    float jValue3 = InterpolateCubic(
                        data[j3, i1], data[j3, i2], data[j3, i3], data[j3, i4], iFraction);
                    float jValue4 = InterpolateCubic(
                        data[j4, i1], data[j4, i2], data[j4, i3], data[j4, i4], iFraction);
                    ret[j, i] = InterpolateCubic(
                        jValue1, jValue2, jValue3, jValue4, jFraction);
                }
            }
        });

        return ret;
    }
}