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)
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)
I'm a bit confused on the third degree polynomial used.
Yes it gives the correct values in 0 and 1, but the derivates of neighbouring cells does not fit, as far as I can calculate. If the grid-data is linear, it does not even return a line....
And it is not point symmetric in x=0.5
The polynomial that fits in 0 and 1 AND also have the same derivates for neighbouring cells, and thus is smooth is (almost) as easy to calculate.
(and it reduces to linear form if that fits the data)
//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;
}
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):
static float[,] GetBicubicInterpolation(float[,] data, int outWidth, int outHeight)
{
// interpolate to a value "fraction" between "v1" and "v2"
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;
}
var width = data.GetLength(1);
var height = data.GetLength(0);
var ret = new float[outHeight, outWidth];
for (int j = 0; j < outHeight; ++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;
}
This should be parallelizable...