我浏览了互联网，在双三次插值方面，我找不到一个简单的方程。维基百科关于这个主题的页面没有太大帮助，所以有什么简单的方法可以学习双三次插补的工作原理和实现方法吗？我用它来生成Perlin Noise，但使用双线性插值是满足我需求的一种方法(我已经尝试过了)。
如果有人能通过一个好的网站或只是一个答案为我指明正确的方向，我将不胜感激
使用这个(感谢Ahmet Kakıcı发现了这个)，我找到了如何添加Bicubic插值。对于那些也在寻找答案的人，这里是我使用的：
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插值，我首先得到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 );
其中数据塔定义为：
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)
接受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，以增加一点踢法。
注意，如果数据在一条线上，这个建议会很好地减少到零
注意：上面Nicholas给出的Bicubic公式(作为答案)中有一个错误。通过对正弦曲线进行插值，我发现它是不正确的。
正确的公式是：
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