Бикубическая интерполяция?
Я просмотрел интернет, и с точки зрения бикубической интерполяции, я не могу найти простое уравнение для этого. Страница Википедии на эту тему оказалась не очень полезной, так есть ли простой способ узнать, как работает бикубическая интерполяция и как ее реализовать? Я использую его для генерации Perlin Noise, но использование билинейной интерполяции - способ нестабильный для моих нужд (я уже пробовал).
Если кто-то может указать мне правильное направление либо с помощью хорошего сайта, либо просто с помощью ответа, я был бы очень признателен. (Я использую C# кстати)
2 ответа
Используя это (спасибо Ахмету Какичи, который нашел это), я понял, как добавить бикубическую интерполяцию. Для тех, кто также ищет ответ, вот что я использовал:
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 );
Где 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соответственно)
Взял ответ 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://stackru.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;
}
}
Примечание. Бикубическая формула, данная Николасом выше (как ответ), содержит ошибку. интерполировав синусоиду, я смог увидеть, что это неверно.
Правильная формула:
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.
Я немного запутался в использовании полинома третьей степени.
Да, это дает правильные значения в 0 и 1, но производные соседних ячеек не подходят, насколько я могу рассчитать. Если данные сетки являются линейными, они даже не возвращают строку....
И это не точка симметрично в х =0,5
Полином, который соответствует 0 и 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 continous smooth curve.
//And we also wants the it to reduce nicely to a line, if that matches the data
//P(x)=Ax^3+Bx^2+Cx-D=((Ax+B)x+C)X+D
//P(0)=D =v0
//P(1)=A+B+C+D =Vp1
//P'(0)=C =(vp1-vm1)/2
//p'(1)=3A+2B+C=(vp2-v0 )/2
//Subtracting first and third from the second
//A+B =vp1-C-D = (vp1+vm1)/2 - v0
//Subtracting that twice and a C from the last
//A=(vp2-v0)/2 - 2(A+B) -C =(vp2-v0)/2 - (Vp1+vm1-2v0) - (vp1-vm1)/2
// = 3(v0-vp1)/2 + (vp2-vm1)/2
//B=(A+B)-A = (vp1+vm1)/2 - v0 - (3(v0-vp1)/2 + (vp2-vm1)/2)
// = vm1 + 2vp1 - (5v0+vp2)/2;
double C = (vp1 - vm1) / 2;
double ApB =vp1 -C -v0;
double A = (vp2 - v0) / 2 - 2 * ApB - C;
double B = ApB - A;
//double B = vm1 + 2 * vp1 - (5 * v0 + vp2) / 2;
//double A = (3*(v0 - vp1) + (vp2 - vm1)) / 2;
return ((A * frac + B) * frac + C) * frac + v0;
}