c# - 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

使用这个（感谢发现这个的 Ahmet Kakıcı），我想出了如何添加双三次插值。对于那些也在寻找答案的人，这是我使用的：

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

为了得到 2D Interpolation，我首先得到 x，然后对 y 进行插值。例如。

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

其中 ndata 定义为：

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 和 inty 是请求坐标的底值。fracx 和 fracy 是输入坐标的小数部分，分别是x-intx和y-inty）

Answer 2

接受 Eske Rahn 的回答并打了一个电话（请注意，下面的代码使用 (j, i) 的矩阵尺寸约定而不是 (x, y) 的图像，但这对于插值而言无关紧要）：

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

Answer 3

我对使用的三次多项式有点困惑。

是的，它在 0 和 1 中给出了正确的值，但据我计算，相邻单元格的导数不适合。如果网格数据是线性的，它甚至不会返回一条线......

它在 x=0.5 中不是点对称的

适合 0 和 1 AND 的多项式对于相邻单元格也具有相同的导数，因此是平滑的，（几乎）同样易于计算。

（如果适合数据，它会简化为线性形式）

在标签中 p 和 m 是加号和减号的简写，例如 vm1 是 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;
    }

受下方ManuTOO评论的启发，我也将其设为四阶，带有可选参数，您可以随意设置/计算，而不会破坏曲线的平滑度。在计算中一路加项基本相同。如果您将 E 设置为零，它将与上述相同。（显然，如果数据实际上在一条线上，则您对 E 的计算必须确保 E 为零才能获得线性输出）

        //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：引用以下@ManTOO 的建议：

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

m_BicubicSharpness 为 1.5，非常接近 Photoshop 的 Bicubic Sharper ；就个人而言，我将其设置为 1.75 以获得额外的效果。

请注意，如果数据在一行上，则此建议很好地减少到零

Answer 4

注意：上面尼古拉斯给出的双三次公式（作为答案）有一个错误。通过插值正弦曲线，我发现它不正确。

正确的公式是：

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;

有关推导，请参阅https://www.paulinternet.nl/?page=bicubic