转载自https://lianera.github.io/post/2018/fast-trilinear-interpolation/
快速三线性插值
最近需要对一个体素进行插值,并且应用到一张大图像上。这个本来用三线性插值很容易就实现了,但是体素的尺寸很小,长宽高大概20x15x10的大小,图像的尺寸非常大,差不多是4000x3000,等我实现完普通的三线性插值后,发现算法速度实在太慢,单是这一项三线性插值耗费的时间就占了总程序的一半多。其实如果能够用GPU,利用OpenGL里面的3d texture是最方便,也最快的,但是图像从内存传送到GPU需要耗费一定时间,并且很多情况下没法用GPU。所以只能推导一下,在体素上事先就进行预计算,然后再用参数应用到大图像上。
原理
借用一张常用的三线性插值示意图
三线性插值实际上可以通过7次线性插值实现。首先假定使用左手坐标系,z轴朝上,在x方向上对四条棱边进行四次插值,得到每条边上的四个点的值c00,c01,c10,c11,然后在y方向上对四个点进行插值,得到两条线段以及中间的两个点的值c0,c1,然后在z方向上插值,得到最终的c点的值。
我们知道,对于线性插值来说,其公式为:
其中包括四次加减运算,三次乘除运算,因此总的来说,进行一次三线性插值,运算总数为28次加减运算,21次乘除运算,并且往体素中寻找位置并且取值还要花费大量的运算,对于大图像来说,总的运算耗费的时间非常大。
假定立方体原点为(x0,y0,z0),离原点最远的点为(x1,y1,z1),立方体的每个顶点的值在体素中是知道的。值得注意的是,虽然现在计算的坐标位置实际上是插值时的位置,但是实际上插值只与待求点(x,y,z)与(x0,y0,z0)和(x1,y1,z1)的相对位置比有关,不过为了后续插值的方便,建议首先将点映射到目标坐标上,这样在插值的时候就不需要再进行坐标映射。首先对x方向上插值,得到yz平面上四个点的值分别为f1,f2,f3,f4:
注意一下,这里面隐含了x1和x0不能相等,如果两者相等,则f1=f(x0,y0,z0)=f(x1,y0,z0)恒成立,所以直接设置成x1=1,x0=0,再用上述公式即可。注意到,这里面的矩阵是不包含未知数的,这是个常数矩阵。这样一来,我们可以事先在体素上计算好矩阵,最后再应用到图像上就好了。当然我们这还只是一次线性插值,现在继续插值。我们令上面式子f1,f2,f3,f4中右侧矩阵分别为M1,M2,M3,M4,因为是常数矩阵,我们将他们用常数a,b替换。同样利用线性插值方法,现在在y轴方向上插值,计算1,2插值的点5,以及3,4插值的点6
实现
预计算参数的实现直接按照上面的常数参数计算即可,因为预计算是在小的体素上面实现的,所以不用太担心速度。
c++
/* precompute trilinear coefficients
* v: f(x0,y0,z0), f(x1,y0,z0), f(x0,y1,z0), f(x1,y1,z0)
* f(x0,y0,z1), f(x1,y0,z1), f(x0,y1,z1), f(x1,y1,z1)
*/
template <typename T>
array<T,8> PreComp(float x0, float y0, float z0,
float x1, float y1, float z1, const array<T,8>& v)
{
const float epsilon = 10e-6f;
if (x1 - x0 < epsilon) { x1 = 1.f; x0 = 0.f; }
if (y1 - y0 < epsilon) { y1 = 1.f; y0 = 0.f; }
if (z1 - z0 < epsilon) { z1 = 1.f; z0 = 0.f; }
float deno = (x1 - x0)*(y1 - y0)*(z1 - z0);
float nume = 1.f / deno;
T a1 = x1 * v[0] - x0 * v[1];
T a2 = x1 * v[2] - x0 * v[3];
T a3 = x1 * v[4] - x0 * v[5];
T a4 = x1 * v[6] - x0 * v[7];
T b1 = v[1] - v[0];
T b2 = v[3] - v[2];
T b3 = v[5] - v[4];
T b4 = v[7] - v[6];
T c5 = y1 * a1 - y0 * a2;
T c6 = y1 * a3 - y0 * a4;
T d5 = a2 - a1;
T d6 = a4 - a3;
T e5 = y1 * b1 - y0 * b2;
T e6 = y1 * b3 - y0 * b4;
T f5 = b2 - b1;
T f6 = b4 - b3;
T h1 = z1 * c5 - z0 * c6;
T h2 = z1 * d5 - z0 * d6;
T h3 = z1 * e5 - z0 * e6;
T h4 = z1 * f5 - z0 * f6;
T h5 = c6 - c5;
T h6 = d6 - d5;
T h7 = e6 - e5;
T h8 = f6 - f5;
array<T, 8> h = {h1*nume, h2*nume, h3*nume, h4*nume,
h5*nume, h6*nume, h7*nume, h8*nume};
return h;
}
而插值的时候,速度就是比较重要的了,因为大体素与小体素对应的时候,block是共享参数的,并且即使是变量x,y,z也可以预先计算一部分,不用每个cell都计算。具体来说,插值的时候直接应用f7的计算公式即可:
T val = zf * (yf*(h[7] * xf + h[5]) + h[6] * xf + h[4])
+ yf * (h[3] * xf + h[1]) + h[2] * xf + h[0];
结果
对一个20x20x20的小体素进行预计算,然后用上述的三线性插值放大5倍,形成一个100x100x100大小的体素。
原体素
完整代码:https://lianera.github.io/post/2018/fast-trilinear-interpolation/main.cpp