programa - Interpolación bilineal en C/C++ y CUDA

metodo de lagrange codigo c++

Question

Quiero emular el comportamiento de la interpolación bilineal de CUDA en la CPU, pero descubrí que el valor de retorno de tex2D no parece ajustarse a la fórmula bilineal .

Supongo que la conversión de los coeficientes de interpolación del formato de coma float al formato de punto fijo de 8 bits con 8 bits de valor fraccionario [1] da como resultado valores diferentes.

De acuerdo con la fórmula de conversión [2, línea 106] , el resultado de la conversión será el mismo que el float entrada cuando el coeficiente es 1/2 1/2^n , con n=0,1,..., 8 , pero yo todavía (no siempre) recibe valores extraños.

A continuación, presento un ejemplo de valores extraños. En este caso, siempre ocurren valores raros cuando id = 2*n+1 , ¿alguien podría decirme por qué?

Src Array:

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

Definición de Textura:

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

Función Kernel:

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;

}

Resultado Bilineal de 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

Resultado de la 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

Dejo un código simple en mi github [3] , después de ejecutar el programa obtendrás dos archivos en D:/ .

Editar 2014/01/20

tex2D el programa con diferentes incrementos y encuentro la especificación de tex2D "cuando alpha multiplicada beta es menor que 0.00390625 , el retorno de tex2D no coincide con la fórmula de interpolación bilineal"

Answer 1

La fórmula incorrecta de interpolación bilineal hace que la textura sea extraña.

Fórmula - 1: puedes encontrarla en el apéndice de cuda o wiki fácilmente

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

Fórmula - 2: reducir los tiempos de multiplicar

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])

Si utiliza el formato de punto fijo de 9 bits para la fórmula 1, obtendrá el resultado de la coincidencia de la búsqueda de textura, pero la fórmula 2 funciona bien.

Conclusión
Si desea emular la interpolación bilineal implementada por cuda texture, debe usar la fórmula 3. ¡Pruébelo!

Fórmula - 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;
}

Answer 2

Los interpolantes UV se truncan a 9 bits, no los valores de téxel participantes. En el Capítulo 10 (Texturizado) del Manual de CUDA, esto se describe en detalle (incluido el código de emulación de CPU) para el caso de 1D. El código es de código abierto y se puede encontrar en https://github.com/ArchaeaSoftware/cudahandbook/blob/master/texturing/tex1d_9bit.cu

Answer 3

Ya se han proporcionado respuestas satisfactorias a esta pregunta, por lo que ahora solo quiero ofrecer un compendio de información útil sobre la interpolación bilineal, cómo se puede implementar en C ++ y las diferentes formas en que se puede hacer en CUDA.

Matemáticas detrás de la interpolación bilineal

Supongamos que la función original T(x, y) se muestrea en la cuadrícula regular cartesiana de puntos (i, j) con 0 <= i < M1 , 0 <= j < M2 e i y j enteros. Para cada valor de y , uno puede usar primero 0 <= a < 1 para representar un punto arbitrario i + a comprendido entre i e i + 1 . Luego, se puede realizar una interpolación lineal a lo largo del eje y = j (que es paralela al eje x ) en ese punto

donde r(x,y) es la función que interpola las muestras de T(x,y) . Lo mismo se puede hacer para la línea y = j + 1 , obteniendo

Ahora, para cada i + a , se puede realizar una interpolación a lo largo del eje y en las muestras r(i+a,j) r(i+a,j+1) . En consecuencia, si uno usa 0 <= b < 1 para representar un punto arbitrario j + b ubicado entre j y j + 1 , entonces una interpolación lineal a lo largo del eje x = i + a (que es paralela al eje y ) puede ser funcionó, por lo que obtener el resultado final

Tenga en cuenta que las relaciones entre i , j , a , b , y son las siguientes

Implementación C / C ++

Permítanme enfatizar que esta implementación, así como las siguientes CUDA, asumen, como se hizo al principio, que las muestras de T están ubicadas en la grilla cartesiana de puntos (i, j) con 0 <= i < M1 , 0 <= j < M2 e i y j enteros (espaciado entre unidades). Además, la rutina se proporciona en aritmética compleja ( float2 ) de precisión simple, pero se puede moldear fácilmente en otras aritméticas de interés.

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);

        }   
    }
}

Las sentencias if/else dentro del código anterior son simplemente comprobaciones de límites. Si la muestra cae fuera de [0, M1-1] x [0, M2-1] , entonces se establece en 0 .

Implementación estándar de CUDA

Esta es una implementación CUDA "estándar" que rastrea la CPU uno anterior. Sin uso de memoria de textura.

__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);

   } 
}

Implementación CUDA con búsqueda de textura

Esta es la misma implementación anterior, pero ahora se accede a la memoria global mediante la memoria caché de texturas. Por ejemplo, se accede a T[i,j] como

tex2D(d_texture_fetch_float,ind_x,ind_y);

(donde, por supuesto, ind_x = i e ind_y = j , y d_texture_fetch_float se supone que es una variable de ámbito global) en lugar de

d_data[ind_y*M1+ind_x];

Tenga en cuenta que las capacidades de filtrado de texturas cableadas no se explotan aquí. La rutina siguiente tiene la misma precisión que la anterior y podría resultar algo más rápida que la de las arquitecturas CUDA antiguas.

__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);

   } 
}

La unión de textura se puede hacer de acuerdo a

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));

}

Tenga en cuenta que ahora no necesitamos ninguna verificación de límites if/else , porque la textura fijará automáticamente a cero las muestras que caen fuera de la región de muestreo [0, M1-1] x [0, M2-1] , gracias a las instrucciones

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

Implementación de CUDA con interpolación de texturas

Esta es la última implementación y utiliza las capacidades de filtrado de texturas.

__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); }
}

Tenga en cuenta que la fórmula de interpolación implementada por esta función es la misma que la derivada anteriormente, pero ahora

donde x_B = x - 0.5 y y_B = y - 0.5 . Esto explica el desplazamiento 0.5 en la instrucción

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

En este caso, la unión de texturas debe hacerse de la siguiente manera

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));

}

Tenga en cuenta que, como ya se mencionó en las otras respuestas, a y b se almacenan en formato de punto fijo de 8 bits con 8 bits de valor fraccionario, por lo que este enfoque será muy rápido, pero menos preciso que los anteriores.