Bicubic Interpolation?

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)

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;
    }
}