我浏览了互联网，就双三次插值而言，我找不到一个简单的方程。维基百科关于这个主题的页面不是很有帮助，那么有没有简单的方法来学习双三次插值的工作原理以及如何实现它？我正在使用它来生成 Perlin 噪声，但使用双线性插值可以满足我的需要（我已经尝试过）。
如果有人可以通过一个好的网站或只是一个答案为我指明正确的方向，我将不胜感激。 （顺便说一下，我用的是C#）
【问题讨论】：
使用this（感谢发现这个的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）
【讨论】：
A*Mathf.Pow(fracy,3)+B*Mathf.Pow(fracy,2)+C*fracy+D; 在时间和精度方面效率不高。我建议你这样实现：D + fracy * (C + fracy * (B + fracy * A))
接听 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;
}
}
【讨论】：
我对使用的三次多项式有点困惑。
是的，它在 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;
添加： 以下引用@ManTOO 的建议：
double E = (v0 - vm1 + vp1 - vp2) * m_BicubicSharpness;
m_BicubicSharpness 为 1.5，非常接近 Photoshop 的 Bicubic Sharper ；就个人而言，我将其设置为 1.75 以获得额外的效果。
请注意，如果数据在一行上，则此建议会很好地减少到零
【讨论】：
注意：上面尼古拉斯给出的双三次公式（作为答案）有一个错误。通过插值正弦曲线，我发现它不正确。
正确的公式是：
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;
【讨论】：