Cubic interpolation立方插值
来源:互联网 发布:不属于c语言标识符 编辑:程序博客网 时间:2024/05/16 05:09
原文地址:http://www.paulinternet.nl/?page=bicubic
On this page you can find explanation about (n-)cubic interpolation and implementations in Java and C++.Anything at this page may be copied and modified.
Please contact me if you find an error.
Cubic interpolation
If the values of a function f(x) and its derivative are known at x=0 and x=1, then the function can be interpolated on the interval [0,1] using a third degree polynomial. This is called cubic interpolation. The formula of this polynomial can be easily derived.
A third degree polynomial and its derivative:
For the green curve:
The values of the polynomial and its derivative at x=0 and x=1:
The four equations above can be rewritten to this:
And there we have our cubic interpolation formula.
Interpolation is often used to interpolate between a list of values. In that case we don't know the derivative of the function. We could simply use derivative 0 at every point, but we obtain smoother curves when we use the slope of a line between the previous and the next point as the derivative at a point. In that case the resulting polynomial is called a Catmull-Rom spline. Suppose you have the values p0, p1, p2 and p3 at respectively x=-1, x=0, x=1, and x=2. Then we can assign the values of f(0), f(1), f'(0) and f'(1) using the formulas below to interpolate between p1 and p2.
Combining the last four formulas and the preceding four, we get:
So our cubic interpolation formula becomes:
The first and the last interval
We used the two points left of the interval and the two points right of the inverval as inputs for the interpolation function. But what if we want to interpolate between the first two or last two elements of a list? Then we have no p0 or no p3. The solution is to imagine an extra point at each end of the list. In other words, we have to make up a value for p0 and p3 when interpolating the leftmost and rightmost interval respectively. Two ways to do this are:
- Repeat the first and the last point.
Left: p0 = p1
Right: p3 = p2 - Let the end point be in the middle of a line between the imaginary point and the point next to the end point.
Left: p0 = 2p1 - p2
Right: p3 = 2p2 - p1
Bicubic interpolation
Bicubic interpolation is cubic interpolation in two dimensions. I'll only consider the case where we want to interpolate a two dimensional grid. We can use the cubic interpolation formula to construct the bicubic interpolation formula.
Suppose we have the 16 points pij, with i and j going from 0 to 3 and with pij located at (i-1, j-1). Then we can interpolate the area [0,1] x [0,1] by first interpolating the four columns and then interpolating the results in the horizontal direction. The formula becomes:
Bicubic interpolation can be used to resize images. However, this is (currently) out of the scope of this article. Please don't ask me questions about that.
Java implementation
public class CubicInterpolator{public static double getValue (double[] p, double x) {return p[1] + 0.5 * x*(p[2] - p[0] + x*(2.0*p[0] - 5.0*p[1] + 4.0*p[2] - p[3] + x*(3.0*(p[1] - p[2]) + p[3] - p[0])));}}public class BicubicInterpolator extends CubicInterpolator{private double[] arr = new double[4];public double getValue (double[][] p, double x, double y) {arr[0] = getValue(p[0], y);arr[1] = getValue(p[1], y);arr[2] = getValue(p[2], y);arr[3] = getValue(p[3], y);return getValue(arr, x);}}public class TricubicInterpolator extends BicubicInterpolator{private double[] arr = new double[4];public double getValue (double[][][] p, double x, double y, double z) {arr[0] = getValue(p[0], y, z);arr[1] = getValue(p[1], y, z);arr[2] = getValue(p[2], y, z);arr[3] = getValue(p[3], y, z);return getValue(arr, x);}}
C++ implementation
#include <iostream>#include <assert.h>double cubicInterpolate (double p[4], double x) {return p[1] + 0.5 * x*(p[2] - p[0] + x*(2.0*p[0] - 5.0*p[1] + 4.0*p[2] - p[3] + x*(3.0*(p[1] - p[2]) + p[3] - p[0])));}double bicubicInterpolate (double p[4][4], double x, double y) {double arr[4];arr[0] = cubicInterpolate(p[0], y);arr[1] = cubicInterpolate(p[1], y);arr[2] = cubicInterpolate(p[2], y);arr[3] = cubicInterpolate(p[3], y);return cubicInterpolate(arr, x);}double tricubicInterpolate (double p[4][4][4], double x, double y, double z) {double arr[4];arr[0] = bicubicInterpolate(p[0], y, z);arr[1] = bicubicInterpolate(p[1], y, z);arr[2] = bicubicInterpolate(p[2], y, z);arr[3] = bicubicInterpolate(p[3], y, z);return cubicInterpolate(arr, x);}double nCubicInterpolate (int n, double* p, double coordinates[]) {assert(n > 0);if (n == 1) {return cubicInterpolate(p, *coordinates);}else {double arr[4];int skip = 1 << (n - 1) * 2;arr[0] = nCubicInterpolate(n - 1, p, coordinates + 1);arr[1] = nCubicInterpolate(n - 1, p + skip, coordinates + 1);arr[2] = nCubicInterpolate(n - 1, p + 2*skip, coordinates + 1);arr[3] = nCubicInterpolate(n - 1, p + 3*skip, coordinates + 1);return cubicInterpolate(arr, *coordinates);}}int main () {// Create arraydouble p[4][4] = {{1,3,3,4}, {7,2,3,4}, {1,6,3,6}, {2,5,7,2}};// Interpolatestd::cout << bicubicInterpolate(p, 0.1, 0.2) << '\n';// Or use the nCubicInterpolate functiondouble co[2] = {0.1, 0.2};std::cout << nCubicInterpolate(2, (double*) p, co) << '\n';}
Bicubic interpolation polynomial
This section provides an alternative way to calculate bicubic interpolation. For most purposes this way is probably less practical and efficient than the way it is done above.
We can rewrite the formula for bicubic interpolation as a multivariate polynomial:
With these values for aij, the coefficients:
In Java code we can write this as:
public class CachedBicubicInterpolator{private double a00, a01, a02, a03;private double a10, a11, a12, a13;private double a20, a21, a22, a23;private double a30, a31, a32, a33;public void updateCoefficients (double[][] p) {a00 = p[1][1];a01 = -.5*p[1][0] + .5*p[1][2];a02 = p[1][0] - 2.5*p[1][1] + 2*p[1][2] - .5*p[1][3];a03 = -.5*p[1][0] + 1.5*p[1][1] - 1.5*p[1][2] + .5*p[1][3];a10 = -.5*p[0][1] + .5*p[2][1];a11 = .25*p[0][0] - .25*p[0][2] - .25*p[2][0] + .25*p[2][2];a12 = -.5*p[0][0] + 1.25*p[0][1] - p[0][2] + .25*p[0][3] + .5*p[2][0] - 1.25*p[2][1] + p[2][2] - .25*p[2][3];a13 = .25*p[0][0] - .75*p[0][1] + .75*p[0][2] - .25*p[0][3] - .25*p[2][0] + .75*p[2][1] - .75*p[2][2] + .25*p[2][3];a20 = p[0][1] - 2.5*p[1][1] + 2*p[2][1] - .5*p[3][1];a21 = -.5*p[0][0] + .5*p[0][2] + 1.25*p[1][0] - 1.25*p[1][2] - p[2][0] + p[2][2] + .25*p[3][0] - .25*p[3][2];a22 = p[0][0] - 2.5*p[0][1] + 2*p[0][2] - .5*p[0][3] - 2.5*p[1][0] + 6.25*p[1][1] - 5*p[1][2] + 1.25*p[1][3] + 2*p[2][0] - 5*p[2][1] + 4*p[2][2] - p[2][3] - .5*p[3][0] + 1.25*p[3][1] - p[3][2] + .25*p[3][3];a23 = -.5*p[0][0] + 1.5*p[0][1] - 1.5*p[0][2] + .5*p[0][3] + 1.25*p[1][0] - 3.75*p[1][1] + 3.75*p[1][2] - 1.25*p[1][3] - p[2][0] + 3*p[2][1] - 3*p[2][2] + p[2][3] + .25*p[3][0] - .75*p[3][1] + .75*p[3][2] - .25*p[3][3];a30 = -.5*p[0][1] + 1.5*p[1][1] - 1.5*p[2][1] + .5*p[3][1];a31 = .25*p[0][0] - .25*p[0][2] - .75*p[1][0] + .75*p[1][2] + .75*p[2][0] - .75*p[2][2] - .25*p[3][0] + .25*p[3][2];a32 = -.5*p[0][0] + 1.25*p[0][1] - p[0][2] + .25*p[0][3] + 1.5*p[1][0] - 3.75*p[1][1] + 3*p[1][2] - .75*p[1][3] - 1.5*p[2][0] + 3.75*p[2][1] - 3*p[2][2] + .75*p[2][3] + .5*p[3][0] - 1.25*p[3][1] + p[3][2] - .25*p[3][3];a33 = .25*p[0][0] - .75*p[0][1] + .75*p[0][2] - .25*p[0][3] - .75*p[1][0] + 2.25*p[1][1] - 2.25*p[1][2] + .75*p[1][3] + .75*p[2][0] - 2.25*p[2][1] + 2.25*p[2][2] - .75*p[2][3] - .25*p[3][0] + .75*p[3][1] - .75*p[3][2] + .25*p[3][3];}public double getValue (double x, double y) {double x2 = x * x;double x3 = x2 * x;double y2 = y * y;double y3 = y2 * y;return (a00 + a01 * y + a02 * y2 + a03 * y3) + (a10 + a11 * y + a12 * y2 + a13 * y3) * x + (a20 + a21 * y + a22 * y2 + a23 * y3) * x2 + (a30 + a31 * y + a32 * y2 + a33 * y3) * x3;}}
This implementation can run faster than the implementation above when getValue is called multiple times for one call to updateCoefficients.
Similar implementations could be written for other dimensions than two.
- Cubic interpolation立方插值
- cubic convolution interpolation (三次卷积插值)
- 插值(Interpolation)
- 插值(interpolation)
- 插值(Interpolation)
- stylus之插值(Interpolation)
- 立方插值
- 立方插值
- 重心座标插值(Barycentric Interpolation)
- OpenGL 实现Interpolation插值算法
- Interpolation Search 插值搜索算法
- less-Variable Interpolation(变量插值)
- scipy module 1 -- Scipy Tutorial 2 Interpolation(插值)
- 视频插值--Video Frame Interpolation via Adaptive Separable Convolution
- 图像放缩------双立方插值
- 图片双立方插值C++版本
- 图像缩放双立方插值
- Opencv 三次样条曲线(Cubic Spline)插值
- 非标准的WINDOW API 根据路径遍历文件 c++
- 【精卫四轴】开源制作过程
- GRE写作必备句型
- 编写一个类,完成Round类实现的功能。(圆的面积,周长,圆柱体积)
- Spring -- ApplicationContext使用及Bean的配置
- Cubic interpolation立方插值
- day05--面向对象--String类
- 2012年5月SAT香港真题解析
- POJ2429 GCD & LCM Inverse pollard_rho大整数分解
- HttpServletRequest对象方法的用法
- 正则表达式怎么匹配非汉字非字母非数字之外的字符
- OC05_NSString(15/7/20)
- iOS事件机制,以及不同手势使用touchesBegan等表现形式
- Python变量的作用域