C/C++和CUDA中的双线性插值

Question

我想在CPU上模拟CUDA双线性插值的行为，但我发现tex2D的返回值似乎不符合双线性公式。

我猜测将插值系数从float转换为具有8位小数位的9位定点格式[1]会导致不同的值。

根据转换公式[2, 第106行]，当系数为1/2^n时，n=0,1,..., 8，转换结果与输入float相同，但我仍然（不总是）收到奇怪的值。

下面我报告一个奇怪值的例子。 在这种情况下，当id = 2 * n + 1时总是出现奇怪的值，请问有人能告诉我原因吗？

源数组：

Src[0][0] =  38;  
Src[1][0] =  39;  
Src[0][1] = 118;  
Src[1][1] =  13;  

纹理定义：

static texture<float4, 2, cudaReadModeElementType> texElnt;
texElnt.addressMode[0] = cudaAddressModeClamp;
texElnt.addressMode[1] = cudaAddressModeClamp;
texElnt.filterMode = cudaFilterModeLinear;
texElnt.normalized = false;

核函数：

static __global__ void kernel_texElnt(float* pdata, int w, int h, int c, float stride/*0.03125f*/) {
    const int gx = blockIdx.x*blockDim.x + threadIdx.x;
    const int gy = blockIdx.y*blockDim.y + threadIdx.y;
    const int gw = gridDim.x * blockDim.x;
    const int gid = gy*gw + gx;
    if (gx >= w || gy >= h) {
        return;
    }

    float2 pnt;
    pnt.x = (gx)*(stride)/*1/32*/;
    pnt.y = 0.0625f/*1/16*/;

    float4 result = tex2D( texElnt, pnt.x + 0.5, pnt.y + 0.5f);
    pdata[gid*3 + 0] = pnt.x;
    pdata[gid*3 + 1] = pnt.y;
    pdata[gid*3 + 2] = result.x;

}

CUDA的双线性计算结果

id  pnt.x   pnt.y   tex2D
 0  0.00000 0.0625  43.0000000  
 1  0.03125 0.0625  42.6171875  
 2  0.06250 0.0625  42.6484375  
 3  0.09375 0.0625  42.2656250  
 4  0.12500 0.0625  42.2968750  
 5  0.15625 0.0625  41.9140625  
 6  0.18750 0.0625  41.9453125  
 7  0.21875 0.0625  41.5625000  
 8  0.25000 0.0625  41.5937500  
 9  0.28125 0.0625  41.2109375  
 0  0.31250 0.0625  41.2421875  
10  0.34375 0.0625  40.8593750  
11  0.37500 0.0625  40.8906250  
12  0.40625 0.0625  40.5078125  
13  0.43750 0.0625  40.5390625  
14  0.46875 0.0625  40.1562500  
15  0.50000 0.0625  40.1875000  
16  0.53125 0.0625  39.8046875  
17  0.56250 0.0625  39.8359375  
18  0.59375 0.0625  39.4531250  
19  0.62500 0.0625  39.4843750  
20  0.65625 0.0625  39.1015625  
21  0.68750 0.0625  39.1328125  
22  0.71875 0.0625  38.7500000  
23  0.75000 0.0625  38.7812500  
24  0.78125 0.0625  38.3984375  
25  0.81250 0.0625  38.4296875  
26  0.84375 0.0625  38.0468750  
27  0.87500 0.0625  38.0781250  
28  0.90625 0.0625  37.6953125  
29  0.93750 0.0625  37.7265625  
30  0.96875 0.0625  37.3437500  
31  1.00000 0.0625  37.3750000

CPU结果：

// convert coefficient ((1-α)*(1-β)), (α*(1-β)), ((1-α)*β), (α*β) to fixed point format  

id  pnt.x   pnt.y   tex2D
 0  0.00000 0.0625 43.00000000  
 1  0.03125 0.0625 43.23046875  
 2  0.06250 0.0625 42.64843750  
 3  0.09375 0.0625 42.87890625  
 4  0.12500 0.0625 42.29687500  
 5  0.15625 0.0625 42.52734375  
 6  0.18750 0.0625 41.94531250  
 7  0.21875 0.0625 42.17578125  
 8  0.25000 0.0625 41.59375000  
 9  0.28125 0.0625 41.82421875  
 0  0.31250 0.0625 41.24218750  
10  0.34375 0.0625 41.47265625  
11  0.37500 0.0625 40.89062500  
12  0.40625 0.0625 41.12109375  
13  0.43750 0.0625 40.53906250  
14  0.46875 0.0625 40.76953125  
15  0.50000 0.0625 40.18750000  
16  0.53125 0.0625 40.41796875  
17  0.56250 0.0625 39.83593750  
18  0.59375 0.0625 40.06640625  
19  0.62500 0.0625 39.48437500  
20  0.65625 0.0625 39.71484375  
21  0.68750 0.0625 39.13281250  
22  0.71875 0.0625 39.36328125  
23  0.75000 0.0625 38.78125000  
24  0.78125 0.0625 39.01171875  
25  0.81250 0.0625 38.42968750  
26  0.84375 0.0625 38.66015625  
27  0.87500 0.0625 38.07812500  
28  0.90625 0.0625 38.30859375  
29  0.93750 0.0625 37.72656250  
30  0.96875 0.0625 37.95703125  
31  1.00000 0.0625 37.37500000

我在我的 GitHub [3]上放置了一份简单代码。运行程序后，您将在D:\目录下找到两个文件。

编辑于2014/01/20

我使用不同的增量运行了该程序，并发现tex2D的规格为"当alpha乘以beta小于0.00390625时，tex2D的返回值与双线性插值公式不匹配"

- user1995868

---

你能否添加一个最短的完整示例，以便其他人可以编译和运行？ - talonmies

感谢您的建议@talonmies，我提供了一个示例代码的链接。 - user1995868

Answer 1

已经有令人满意的答案回答了这个问题，现在我只想提供一些关于双线性插值的有用信息，包括如何在C++中实现以及在CUDA中不同的实现方式。

双线性插值背后的数学原理

假设原始函数T(x, y)在笛卡尔坐标系的规则网格点(i, j)上进行采样，其中0 <= i < M1，0 <= j < M2且i和j为整数。对于每个y的值，可以使用0 <= a < 1来表示位于i和i + 1之间的任意点i + a。然后，在该点沿着与x轴平行的y = j轴进行线性插值，得到

其中r(x,y)是插值T(x,y)样本的函数。同样，可以在线y = j + 1上执行相同的操作，得到

现在，对于每个i + a，可以在样本r(i+a,j)和r(i+a,j+1)上沿着y轴进行插值。因此，如果使用0 <= b < 1来表示位于j和j + 1之间的任意点j + b，则可以沿着与y轴平行的x = i + a轴进行线性插值，从而得到最终结果

请注意，i、j、a、b、x和y之间的关系如下

C/C++实现

请注意，这个实现以及下面的CUDA实现都假定，就像一开始做的那样，T的样本位于点j的笛卡尔正则网格上，其中的i和j是整数（单位间距），其范围为0 <= i < M1，0 <= j < M2。此外，该程序使用单精度、复杂（float2）算法来提供服务，但可以轻松转换为其他算法。请保留HTML标签。

void bilinear_interpolation_function_CPU(float2 * __restrict__ h_result, float2 * __restrict__ h_data, 
                                         float * __restrict__ h_xout, float * __restrict__ h_yout, 
                                         const int M1, const int M2, const int N1, const int N2){

    float2 result_temp1, result_temp2;
    for(int k=0; k<N2; k++){
        for(int l=0; l<N1; l++){

            const int   ind_x = floor(h_xout[k*N1+l]); 
            const float a     = h_xout[k*N1+l]-ind_x; 

            const int   ind_y = floor(h_yout[k*N1+l]); 
            const float b     = h_yout[k*N1+l]-ind_y; 

            float2 h00, h01, h10, h11;
            if (((ind_x)   < M1)&&((ind_y)   < M2)) h00 = h_data[ind_y*M1+ind_x];       else    h00 = make_float2(0.f, 0.f);
            if (((ind_x+1) < M1)&&((ind_y)   < M2)) h10 = h_data[ind_y*M1+ind_x+1];     else    h10 = make_float2(0.f, 0.f);
            if (((ind_x)   < M1)&&((ind_y+1) < M2)) h01 = h_data[(ind_y+1)*M1+ind_x];   else    h01 = make_float2(0.f, 0.f);
            if (((ind_x+1) < M1)&&((ind_y+1) < M2)) h11 = h_data[(ind_y+1)*M1+ind_x+1]; else    h11 = make_float2(0.f, 0.f);

            result_temp1.x = a * h10.x + (-h00.x * a + h00.x); 
            result_temp1.y = a * h10.y + (-h00.y * a + h00.y);

            result_temp2.x = a * h11.x + (-h01.x * a + h01.x);
            result_temp2.y = a * h11.y + (-h01.y * a + h01.y);

            h_result[k*N1+l].x = b * result_temp2.x + (-result_temp1.x * b + result_temp1.x);
            h_result[k*N1+l].y = b * result_temp2.y + (-result_temp1.y * b + result_temp1.y);

        }   
    }
}

上述代码中的if/else语句只是边界检查。如果样本超出了[0，M1-1] x [0，M2-1]，则将其设置为0。

标准CUDA实现

这是一个“标准”的CUDA实现，跟踪上面的CPU实现。没有使用纹理内存。

__global__ void bilinear_interpolation_kernel_GPU(float2 * __restrict__ d_result, const float2 * __restrict__ d_data, 
                                                  const float * __restrict__ d_xout, const float * __restrict__ d_yout, 
                                                  const int M1, const int M2, const int N1, const int N2)
{
   const int l = threadIdx.x + blockDim.x * blockIdx.x;
   const int k = threadIdx.y + blockDim.y * blockIdx.y;

   if ((l<N1)&&(k<N2)) {

       float2 result_temp1, result_temp2;

       const int    ind_x = floor(d_xout[k*N1+l]); 
       const float  a     = d_xout[k*N1+l]-ind_x; 

       const int    ind_y = floor(d_yout[k*N1+l]); 
       const float  b     = d_yout[k*N1+l]-ind_y; 

       float2 d00, d01, d10, d11;
       if (((ind_x)   < M1)&&((ind_y)   < M2))  d00 = d_data[ind_y*M1+ind_x];       else    d00 = make_float2(0.f, 0.f);
       if (((ind_x+1) < M1)&&((ind_y)   < M2))  d10 = d_data[ind_y*M1+ind_x+1];     else    d10 = make_float2(0.f, 0.f);
       if (((ind_x)   < M1)&&((ind_y+1) < M2))  d01 = d_data[(ind_y+1)*M1+ind_x];   else    d01 = make_float2(0.f, 0.f);
       if (((ind_x+1) < M1)&&((ind_y+1) < M2))  d11 = d_data[(ind_y+1)*M1+ind_x+1]; else    d11 = make_float2(0.f, 0.f);

       result_temp1.x = a * d10.x + (-d00.x * a + d00.x); 
       result_temp1.y = a * d10.y + (-d00.y * a + d00.y);

       result_temp2.x = a * d11.x + (-d01.x * a + d01.x);
       result_temp2.y = a * d11.y + (-d01.y * a + d01.y);

       d_result[k*N1+l].x = b * result_temp2.x + (-result_temp1.x * b + result_temp1.x);
       d_result[k*N1+l].y = b * result_temp2.y + (-result_temp1.y * b + result_temp1.y);

   } 
}

使用纹理获取实现CUDA

这与上面的实现相同，但现在通过纹理缓存访问全局内存。例如，T[i,j] 的访问方式为

tex2D(d_texture_fetch_float,ind_x,ind_y);

（其中，当然假设ind_x = i和ind_y = j，并且d_texture_fetch_float被认为是全局范围变量），而不是

d_data[ind_y*M1+ind_x];

请注意，此处未利用硬连线纹理过滤功能。下面的例程与上述例程具有相同的精度，并且在旧CUDA架构上可能比其更快。

__global__ void bilinear_interpolation_kernel_GPU_texture_fetch(float2 * __restrict__ d_result, 
                                                                const float * __restrict__ d_xout, const float * __restrict__ d_yout, 
                                                                const int M1, const int M2, const int N1, const int N2)
{
   const int l = threadIdx.x + blockDim.x * blockIdx.x;
   const int k = threadIdx.y + blockDim.y * blockIdx.y;

   if ((l<N1)&&(k<N2)) {

       float2 result_temp1, result_temp2;

       const int    ind_x = floor(d_xout[k*N1+l]); 
       const float  a     = d_xout[k*N1+l]-ind_x; 

       const int    ind_y = floor(d_yout[k*N1+l]); 
       const float  b     = d_yout[k*N1+l]-ind_y; 

       const float2 d00  = tex2D(d_texture_fetch_float,ind_x,ind_y);    
       const float2 d10  = tex2D(d_texture_fetch_float,ind_x+1,ind_y);
       const float2 d11  = tex2D(d_texture_fetch_float,ind_x+1,ind_y+1);
       const float2 d01  = tex2D(d_texture_fetch_float,ind_x,ind_y+1);

       result_temp1.x = a * d10.x + (-d00.x * a + d00.x); 
       result_temp1.y = a * d10.y + (-d00.y * a + d00.y);

       result_temp2.x = a * d11.x + (-d01.x * a + d01.x);
       result_temp2.y = a * d11.y + (-d01.y * a + d01.y);

       d_result[k*N1+l].x = b * result_temp2.x + (-result_temp1.x * b + result_temp1.x);
       d_result[k*N1+l].y = b * result_temp2.y + (-result_temp1.y * b + result_temp1.y);

   } 
}

纹理绑定可以根据以下方式进行：

void TextureBindingBilinearFetch(const float2 * __restrict__ data, const int M1, const int M2)
{
    size_t pitch; 
    float* data_d;
    gpuErrchk(cudaMallocPitch((void**)&data_d,&pitch, M1 * sizeof(float2), M2));
    cudaChannelFormatDesc desc = cudaCreateChannelDesc<float2>();
    gpuErrchk(cudaBindTexture2D(0,&d_texture_fetch_float,data_d,&desc,M1,M2,pitch));
    d_texture_fetch_float.addressMode[0] = cudaAddressModeClamp;
    d_texture_fetch_float.addressMode[1] = cudaAddressModeClamp;
    gpuErrchk(cudaMemcpy2D(data_d,pitch,data,sizeof(float2)*M1,sizeof(float2)*M1,M2,cudaMemcpyHostToDevice));

}

请注意，现在我们不需要进行if/else边界检查，因为纹理将自动夹紧到[0, M1-1] x [0, M2-1]采样区域之外的样本将被裁剪为零，这要感谢指令。

    d_texture_fetch_float.addressMode[0] = cudaAddressModeClamp;
    d_texture_fetch_float.addressMode[1] = cudaAddressModeClamp;

CUDA实现纹理插值

这是最后一种实现方式，利用了纹理过滤的硬件能力。

__global__ void bilinear_interpolation_kernel_GPU_texture_interp(float2 * __restrict__ d_result, 
                                                                 const float * __restrict__ d_xout, const float * __restrict__ d_yout, 
                                                                 const int M1, const int M2, const int N1, const int N2)
{
   const int l = threadIdx.x + blockDim.x * blockIdx.x;
   const int k = threadIdx.y + blockDim.y * blockIdx.y;

   if ((l<N1)&&(k<N2)) { d_result[k*N1+l] = tex2D(d_texture_interp_float, d_xout[k*N1+l] + 0.5f, d_yout[k*N1+l] + 0.5f); }
}

请注意，此功能实现的插值公式与上述导出的公式相同，但现在为 其中 x_B = x - 0.5 且 y_B = y - 0.5。这解释了指令中的 0.5 偏移量。

tex2D(d_texture_interp_float, d_xout[k*N1+l] + 0.5f, d_yout[k*N1+l] + 0.5f)

在这种情况下，纹理绑定应该按照以下方式进行。

void TextureBindingBilinearInterp(const float2 * __restrict__ data, const int M1, const int M2)
{
    size_t pitch; 
    float* data_d;
    gpuErrchk(cudaMallocPitch((void**)&data_d,&pitch, M1 * sizeof(float2), M2));
    cudaChannelFormatDesc desc = cudaCreateChannelDesc<float2>();
    gpuErrchk(cudaBindTexture2D(0,&d_texture_interp_float,data_d,&desc,M1,M2,pitch));
    d_texture_interp_float.addressMode[0] = cudaAddressModeClamp;
    d_texture_interp_float.addressMode[1] = cudaAddressModeClamp;
    d_texture_interp_float.filterMode = cudaFilterModeLinear;   // --- Enable linear filtering
    d_texture_interp_float.normalized = false;                  // --- Texture coordinates will NOT be normalized
    gpuErrchk(cudaMemcpy2D(data_d,pitch,data,sizeof(float2)*M1,sizeof(float2)*M1,M2,cudaMemcpyHostToDevice));

}

请注意，正如其他答案中已经提到的那样，a和b以9位定点格式存储，并具有8位小数值，因此这种方法速度非常快，但精度不如上面的方法。

- Vitality

---

如果想进一步探讨如何利用CUDA纹理内存实现双线性插值，请查看编写插值的CUDA核函数。 - Vitality

Answer 2

UV插值器被截断为9位，而不是参与的纹素值。在CUDA手册的第10章（纹理）中，对于1D情况进行了详细描述（包括CPU仿真代码）。代码是开源的，可在https://github.com/ArchaeaSoftware/cudahandbook/blob/master/texturing/tex1d_9bit.cu找到。

- ArchaeaSoftware

---

感谢您的回复，我之前已经阅读了代码（与文章中的链接2相同）。根据链接，如果x = 1/2^n（n = 1, 2，... 8），那么frac应该等于x。但是我使用tex2D得到了不匹配的结果，所以我在这里发布了我的问题。 - user1995868

Answer 3

错误的双线性插值公式会导致纹理获取结果异常。

公式1：您可以在CUDA附录或维基百科中轻松找到它

tex(x,y)=(1−α)(1−β)T[i,j] + α(1−β)T[i+1,j] + (1−α)βT[i,j+1] + αβT[i+1,j+1]

公式 - 2：减少乘法次数

tex(x,y)=T[i,j] + α(T[i+1,j]-T[i,j]) + β(T[i,j+1]-T[i,j]) + αβ(T[i,j]+T[i+1,j+1] - T[i+1, j]-T[i,j+1])

如果您在公式1中使用9位定点格式，会得到纹理获取不匹配的结果，但公式2可以正常工作。

结论： 如果您想模拟由cuda纹理实现的双线性插值，则应使用公式3。试试吧！

公式-3:

tex(x,y)=T[i,j] + frac(α)(T[i+1,j]-T[i,j]) + frac(β)(T[i,j+1]-T[i,j]) + frac(αβ)(T[i,j]+T[i+1,j+1] - T[i+1, j]-T[i,j+1])  

// frac(x) turns float to 9-bit fixed point format with 8 bits of fraction values.     
float frac( float x ) {
    float frac, tmp = x - (float)(int)(x);
    float frac256 = (float)(int)( tmp*256.0f + 0.5f );
    frac = frac256 / 256.0f;
    return frac;
}

- user1995868