BookGPU/Chapters/chapter4/code/levelines_kernels.cu

   1 #define PI 3.1415926
   2
   3 __constant__ float masque[256] ;
   4 /* noyau gaussien 3x3 *//*
   5 __constant__ float noyau[] = {1.0/16,2.0/16,1.0/16,
   6                                                         2.0/16,4.0/16,2.0/16,
   7                                                         1.0/16,2.0/16,1.0/16} ;
   8                                                         __constant__ int rnoyau  = 1 ;*/
   9
  10 /* noyau gaussien 5x5 *//*
  11 __constant__ int noyau[] = {1,4,6,4,1, 4,16,24,16,4, 6,24,36,24,6, 4,16,24,16,4, 1,4,6,4,1} ;
  12 __constant__ int sumNoyau = 256 ;
  13 __constant__ int rnoyau  = 2 ;
  14 __constant__ int NNoyau = 25 ;
  15                                                 */
  16 /* noyau 5x5*/
  17 __constant__ float noyau[] = {1.0/25,1.0/25,1.0/25,1.0/25,1.0/25,
  18                                                           1.0/25,1.0/25,1.0/25,1.0/25,1.0/25,
  19                                                           1.0/25,1.0/25,1.0/25,1.0/25,1.0/25,
  20                                                           1.0/25,1.0/25,1.0/25,1.0/25,1.0/25,
  21                                                           1.0/25,1.0/25,1.0/25,1.0/25,1.0/25} ;
  22 __constant__ int rnoyau  = 2 ;
  23
  24 /* noyau 7x7 *//*
  25 __device__ __constant__ float noyau[] = {1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,
  26                                                                                  1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,
  27                                                                                  1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,
  28                                                                                  1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,
  29                                                                                  1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,
  30                                                                                  1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,
  31                                                                                  1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49} ;
  32 __device__ __constant__ int rnoyau  = 3 ;
  33                            */
  34 //__device__ __constant__ int NNoyau = 49 ;
  35
  36 __device__ __constant__ int noyauH[] = {1,1,1,1,1,1,1} ;
  37 __device__ __constant__ int sumNoyauH = 7 ;
  38 __device__ __constant__ int rnoyauH  = 3 ;
  39 // average 7x7
  40 /*
  41 __device__ __constant__ int noyauV[] = {1,1,1,1,1,1,1} ;
  42 __device__ __constant__ int sumNoyauV = 7 ;
  43 __device__ __constant__ int rnoyauV  = 3 ;
  44 */
  45 // gaussian 5x5
  46 __device__ __constant__ int noyauV[] = {1,4,6,4,1} ;
  47 __device__ __constant__ int rnoyauV  = 2 ;
  48 __constant__ int sumKV = 16 ;
  49 //__device__ __constant__ int noyauV[] = {1,1,1} ;
  50 //__constant__ int sumNoyauV = 16 ;
  51
  52 // declarations des textures
  53 texture<int, 2, cudaReadModeElementType> tex_img_in ;
  54 texture<unsigned short, 2, cudaReadModeElementType> tex_img_ins ;
  55 texture<unsigned char, 2, cudaReadModeElementType> tex_img_inc ;
  56 texture<int, 2, cudaReadModeElementType> tex_img_estim ;
  57 texture<float, 1, cudaReadModeElementType> tex_noyau ;
  58 texture<float, 1, cudaReadModeElementType> tex_noyauConvof ;
  59 texture<short, 1, cudaReadModeElementType> tex_noyauConvos ;
  60 texture<int, 1, cudaReadModeElementType> tex_kern ;
  61
  62
  63 /***************************************************************
  64  * kernel identite pour mesures de debit max effectif possible
  65  ***************************************************************/
  66 __global__ void kernel_ident( unsigned short*output, int i_dim, int j_dim)
  67 {
  68   // coordonnees absolues du point de base
  69   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
  70   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  71   int idx = __mul24(i, j_dim) + j ;              // indice  dans l'image
  72
  73   output[ idx ] = tex2D(tex_img_ins, j, i) ;
  74
  75 }
  76
  77 __global__ void kernel_ident2p( unsigned short  *output, int i_dim, int j_dim)
  78 {
  79   // coordonnees absolues du point de base
  80   int j = __umul24(__umul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
  81   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  82
  83   output[ __umul24(i, j_dim) + j    ] = tex2D(tex_img_inc, j, i) ;
  84   output[ __umul24(i, j_dim) + j +1 ] = tex2D(tex_img_inc, j+1, i) ;
  85
  86 }
  87
  88
  89 __global__ void kernel_ident( unsigned char  *output, int i_dim, int j_dim)
  90 {
  91   // coordonnees absolues du point de base
  92   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
  93   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  94
  95   output[ __umul24(i, j_dim) + j ] = tex2D(tex_img_inc, j, i) ;
  96
  97 }
  98
  99 __global__ void kernel_ident2p( unsigned char  *output, int i_dim, int j_dim)
 100 {
 101   // coordonnees absolues du point de base
 102   int j = __umul24(__umul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
 103   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
 104
 105   output[ __umul24(i, j_dim) + j    ] = tex2D(tex_img_inc, j, i) ;
 106   output[ __umul24(i, j_dim) + j +1 ] = tex2D(tex_img_inc, j+1, i) ;
 107
 108 }
 109
 110 /**
 111  *
 112  * \brief calcule l'estimation initiale
 113  * \author zulu - AND
 114  *
 115  * \param[in] L         Largeur de l'image
 116  * \param[in] i_dim     Hauteur de l'image
 117  * \param[in] j_dim     coté de la fenêtre de moyenneage
 118  *
 119  * \param[out] output   Image estimee initiale
 120  *
 121  * Version texture : l'img originale est supposée en texture.
 122  * L'estimation réalisée correspond a un moyenneur de 'rayon' r
 123  * Execution sur des blocs de threads 2D et une grille 2D
 124  * selon les dimensions de l'image.
 125  *
 126  */
 127
 128 __global__ void kernel_moyenne_shtex( unsigned int *output, unsigned int i_dim, unsigned int j_dim, unsigned int r)
 129 {
 130   int idl, idc ;
 131   //coordonnées ds le bloc
 132   int ib = threadIdx.y ;
 133   int jb = threadIdx.x ;
 134   int idb_row = __mul24(ib,blockDim.x) + jb ;   // indice dans le bloc
 135
 136   // coordonnees absolues du point
 137   int j = __mul24(blockIdx.x,blockDim.x) + jb ;
 138   int i = __mul24(blockIdx.y,blockDim.y) + ib ;
 139   int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
 140   int off = __mul24(r,blockDim.x) ;             // offset du bloc dans le buffer en sh mem
 141
 142
 143   extern __shared__ unsigned int shmem[] ;
 144   volatile unsigned int * buffer = shmem ;
 145   volatile unsigned int * bufferV = &buffer[ blockDim.x*blockDim.y + 2*off ] ;
 146
 147
 148   /***********************************************************************************
 149    *               PASSE HORIZONTALE
 150    ***********************************************************************************/
 151   // charge les pixels de la bande au dessus du bloc
 152   if(ib<r)
 153         {
 154           buffer[ idb_row ] = 0 ;
 155           for (idc=j-r ; idc<=j+r ; idc++ )
 156                 buffer[ idb_row  ] += tex2D(tex_img_in, idc, i-r) ;
 157         }
 158
 159   // charge les pixels de la bande en dessous du bloc
 160   if (ib>=(blockDim.y-r))
 161         {
 162           buffer[ idb_row + 2*off ] = 0 ;
 163           for (idc=j-r ; idc<=j+r ; idc++ )
 164                 buffer[ idb_row + 2*off ] += tex2D(tex_img_in, idc, i+r) ;
 165         }
 166
 167   // charge les pixels du bloc
 168   buffer[ idb_row + off ] = 0 ;
 169   for (idc=j-r ; idc<=j+r ; idc++ )
 170         buffer[ idb_row + off ] += tex2D(tex_img_in, idc, i) ;
 171
 172    __syncthreads() ;
 173   /*****************************************************************************
 174    *                             PASSE VERTICALE
 175    *****************************************************************************/
 176    bufferV[ idb_row ] = 0 ;
 177    for (idl=0 ; idl<=2*r ; idl++ )
 178          bufferV[ idb_row ] += buffer[ (ib+idl)*blockDim.x + jb ] ;
 179
 180   __syncthreads() ;
 181
 182   output[idx] = bufferV[ idb_row ]/((2*r+1)*(2*r+1)) ; //pas juste sur les bords, mais c'est pô grave !
 183
 184 }
 185
 186 /**
 187  *
 188  * \brief calcule la convolution avec un noyau carre
 189  * \author zulu - AND
 190  *
 191  * \param[in] L         Largeur de l'image
 192  * \param[in] i_dim     Hauteur de l'image
 193  * \param[in] j_dim     coté de la fenêtre de moyenneage
 194  *
 195  * \param[out] output   Image estimee initiale
 196  *
 197  * Version texture : l'img originale est supposée en texture.
 198  * Execution sur des blocs de threads 2D et une grille 2D
 199  * correspondant aux dimensions de l'image.
 200  *
 201  */
 202
 203 // convolution non separable
 204 // image en texture et noyau en mem constante
 205 __global__ void kernel_convoGene8( unsigned char  *output, int i_dim, int j_dim)
 206 {
 207   int ic, jc ;
 208   int r = 2 ;
 209   int L=5 ;
 210   float outval0=0.0 ;
 211
 212   // coordonnees absolues du point de base
 213   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
 214   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
 215
 216   for (ic=0 ; ic<L ; ic++)
 217         for(jc=0 ; jc<L ; jc++)
 218           outval0 += masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
 219
 220
 221   output[ __mul24(i, j_dim) + j ] = outval0 ;
 222 }
 223
 224 // convolution non separable 2 pixels/thread
 225 // image en texture et noyau en mem constante
 226 __global__ void kernel_convoGene8x2p( unsigned char  *output, int i_dim, int j_dim )
 227 {
 228   int ic, jc ;
 229   int r = 2 ;
 230   int L=5 ;
 231   float outval0=0.0, outval1=0.0, outval2=0.0 ;
 232
 233   // coordonnees absolues du point de base
 234   int j = __mul24( __mul24(blockIdx.x,blockDim.x) + threadIdx.x, 2) ;
 235   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
 236
 237   for (ic=0 ; ic<L ; ic++)
 238         {
 239           for(jc=1 ; jc<L-1 ; jc++)
 240                 {
 241                   outval0 += masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
 242                 }
 243         }
 244   outval1 = outval0 ;
 245   outval2 = outval0 ;
 246   for (ic=0 ; ic<L ; ic++)
 247         {
 248           outval1 += masque[ __mul24(ic,L) ]*tex2D(tex_img_inc, j-r, i-r+ic) ;
 249           outval2 += masque[ __mul24(ic,L)+L-1 ]*tex2D(tex_img_inc, j-r+L, i-r+ic) ;
 250         }
 251
 252   output[ __mul24(i, j_dim) + j ] = outval1 ;
 253   output[ __mul24(i, j_dim) + j+1 ] = outval2 ;
 254 }
 255
 256 // convolution non separable 2 pixels/thread
 257 // image en texture et noyau en mem constante
 258 __global__ void kernel_convoGene8rx2p( unsigned char  *output, int j_dim, int r )
 259 {
 260   int ic, jc ;
 261   int L=2*r+1 ;
 262   //unsigned char pix ;
 263   float outval0=0.0, outval1=0.0 ;
 264
 265   // coordonnees absolues du point de base
 266   int j = __umul24( __umul24(blockIdx.x,blockDim.x) + threadIdx.x, 2) ;
 267   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
 268   int idx=__umul24(i, j_dim) + j ;
 269
 270   for (ic=0 ; ic<L ; ic++)
 271         {
 272           outval0 += masque[ __umul24(ic,L) ]*tex2D(tex_img_inc, j-r, i-r+ic) ; //bande gauche
 273           //zone commune
 274           for(jc=1 ; jc<L ; jc++)
 275                 {
 276                   //pix = tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
 277                   outval0 += masque[ __umul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
 278                   outval1 += masque[ __umul24(ic,L)+jc-1 ]*tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
 279                 }
 280           outval1 += masque[ __umul24(ic,L)+L-1 ]*tex2D(tex_img_inc, j-r+L, i-r+ic) ; // bande droite
 281         }
 282   output[ idx   ] = outval0 ;
 283   output[ idx+1 ] = outval1 ;
 284 }
 285
 286 // convolution non separable 4 pixels/thread
 287 // image en texture et noyau en mem constante
 288 __global__ void kernel_convoGene8rx4p( unsigned char  *output, int j_dim, int r )
 289 {
 290   int ic, jc ;
 291   int L=2*r+1 ;
 292   unsigned char pix ;
 293   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
 294
 295   // coordonnees absolues du point de base en haut a gauche
 296   int j = __mul24( __umul24( blockIdx.x, blockDim.x) + threadIdx.x, 2) ;
 297   int i = __mul24( __umul24( blockIdx.y, blockDim.y) + threadIdx.y, 2) ;
 298
 299   //zone commune
 300   //forme generique contrib pixel => masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic)
 301   for (ic=1 ; ic<L ; ic++)
 302         {
 303           for(jc=1 ; jc<L ; jc++)
 304                 {
 305                   pix = tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
 306                   outval0 += masque[ __umul24(ic,L)  +jc   ]*pix ;
 307                   outval1 += masque[ __umul24(ic,L)  +jc-1 ]*pix ;
 308                   outval2 += masque[ __umul24(ic-1,L)+jc   ]*pix ;
 309                   outval3 += masque[ __umul24(ic-1,L)+jc-1 ]*pix ;
 310                 }
 311         }
 312   // les blocs peripheriques verticaux
 313   for (ic=1 ; ic<L ; ic++)
 314         {
 315           outval0 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r, i-r+ic) ; // bande A pour b0
 316           outval2 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r, i-r+ic) ; // bande A pour b2
 317
 318           outval1 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+1, i-r+ic) ; //bande C pour b1
 319           outval3 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+1, i-r+ic) ; //bande C pour b3
 320         }
 321   // les blocs peripheriques horizontaux
 322   for (jc=1 ; jc<L ; jc++)
 323         {
 324           outval0 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc, i-r) ; // bande B pour b0
 325           outval1 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc, i-r) ; // bande B pour b1
 326
 327           outval2 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc, i+r+1) ; //bande D pour b2
 328           outval3 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc, i+r+1) ; //bande D pour b3
 329         }
 330
 331   //les coins
 332   {
 333         outval0 += masque[ 0                   ]*tex2D(tex_img_inc, j-r  , i-r) ;   // haut gauche pour  b0
 334         outval1 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+1, i-r) ;   // haut droit pour b1
 335         outval2 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r  , i+r+1) ; // bas gauche pour b2
 336         outval3 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+1, i+r+1) ; //bas droit pour b3
 337   }
 338
 339
 340   // les 2 pixels hauts
 341   output[ __umul24(i, j_dim) + j   ] = outval0 ;
 342   output[ __umul24(i, j_dim) + j+1 ] = outval1 ;
 343   // les 2 pixels bas
 344   output[ __umul24(i+1, j_dim) + j   ] = outval2 ;
 345   output[ __umul24(i+1, j_dim) + j+1 ] = outval3 ;
 346 }
 347
 348 // convolution non separable 4 pixels/thread
 349 // image en texture et noyau en mem constante
 350 __global__ void kernel_convoGene8rx8p( unsigned char  *output, int j_dim, int r )
 351 {
 352   int ic, jc ;
 353   int L=2*r+1 ;
 354   unsigned char pix ;
 355   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
 356   float outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
 357
 358   // coordonnees absolues du point de base en haut a gauche
 359   int j = ( __umul24( blockIdx.x, blockDim.x) + threadIdx.x)<< 2 ;
 360   int i = ( __umul24( blockIdx.y, blockDim.y) + threadIdx.y)<< 1 ;
 361
 362   //zone commune
 363   //forme generique contrib pixel => masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic)
 364   for (ic=1 ; ic<L ; ic++)
 365         {
 366           for(jc=1 ; jc<L ; jc++)
 367                 {
 368                   pix = tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
 369                   outval0 += masque[ __umul24(ic,L)  +jc   ]*pix ;
 370                   outval1 += masque[ __umul24(ic,L)  +jc-1 ]*pix ;
 371                   outval2 += masque[ __umul24(ic-1,L)+jc   ]*pix ;
 372                   outval3 += masque[ __umul24(ic-1,L)+jc-1 ]*pix ;
 373
 374                   pix = tex2D(tex_img_inc, j-r+jc+2, i-r+ic) ;
 375                   outval4 += masque[ __umul24(ic,L)  +jc   ]*pix ;
 376                   outval5 += masque[ __umul24(ic,L)  +jc-1 ]*pix ;
 377                   outval6 += masque[ __umul24(ic-1,L)+jc   ]*pix ;
 378                   outval7 += masque[ __umul24(ic-1,L)+jc-1 ]*pix ;
 379
 380                 }
 381         }
 382   // les blocs peripheriques verticaux
 383   for (ic=1 ; ic<L ; ic++)
 384         {
 385           outval0 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r, i-r+ic) ; // bande A pour b0
 386           outval2 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r, i-r+ic) ; // bande A pour b2
 387
 388           outval1 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+1, i-r+ic) ; //bande C pour b1
 389           outval3 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+1, i-r+ic) ; //bande C pour b3
 390
 391           outval4 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r+2, i-r+ic) ; // bande A pour b0
 392           outval5 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r+2, i-r+ic) ; // bande A pour b2
 393
 394           outval6 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+3, i-r+ic) ; //bande C pour b1
 395           outval7 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+3, i-r+ic) ; //bande C pour b3
 396         }
 397   // les blocs peripheriques horizontaux
 398   for (jc=1 ; jc<L ; jc++)
 399         {
 400           outval0 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc, i-r) ; // bande B pour b0
 401           outval1 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc, i-r) ; // bande B pour b1
 402
 403           outval2 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc, i+r+1) ; //bande D pour b2
 404           outval3 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc, i+r+1) ; //bande D pour b3
 405
 406           outval4 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc+2, i-r) ; // bande B pour b0
 407           outval5 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc+2, i-r) ; // bande B pour b1
 408
 409           outval6 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc+2, i+r+1) ; //bande D pour b2
 410           outval7 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc+2, i+r+1) ; //bande D pour b3
 411         }
 412
 413   //les coins
 414   {
 415         outval0 += masque[ 0                   ]*tex2D(tex_img_inc, j-r  , i-r) ;   // haut gauche pour  b0
 416         outval1 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+1, i-r) ;   // haut droit pour b1
 417         outval2 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r  , i+r+1) ; // bas gauche pour b2
 418         outval3 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+1, i+r+1) ; //bas droit pour b3
 419
 420         outval4 += masque[ 0                   ]*tex2D(tex_img_inc, j-r+2, i-r) ;   // haut gauche pour  b0
 421         outval5 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+3, i-r) ;   // haut droit pour b1
 422         outval6 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r+2, i+r+1) ; // bas gauche pour b2
 423         outval7 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+3, i+r+1) ; //bas droit pour b3
 424   }
 425
 426
 427   // les 2 pixels hauts
 428   output[ __umul24(i, j_dim) + j   ] = outval0 ;
 429   output[ __umul24(i, j_dim) + j+1 ] = outval1 ;
 430   output[ __umul24(i, j_dim) + j+2 ] = outval4 ;
 431   output[ __umul24(i, j_dim) + j+3 ] = outval5 ;
 432   // les 2 pixels bas
 433   output[ __umul24(i+1, j_dim) + j   ] = outval2 ;
 434   output[ __umul24(i+1, j_dim) + j+1 ] = outval3 ;
 435   output[ __umul24(i+1, j_dim) + j+2 ] = outval6 ;
 436   output[ __umul24(i+1, j_dim) + j+3 ] = outval7 ;
 437
 438 }
 439
 440 // convolution non separable 8 pixels/thread
 441 // image en texture et noyau en mem constante
 442 __global__ void kernel_convoGene8x8pL3( unsigned char  *output, int j_dim )
 443 {
 444   int ic, jc ;
 445   int L=3 ;
 446   unsigned char pix ;
 447   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
 448   float outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
 449
 450   // coordonnees absolues du point de base en haut a gauche
 451   int j = ( __umul24( blockIdx.x, blockDim.x) + threadIdx.x)<< 3 ;
 452   int i = ( __umul24( blockIdx.y, blockDim.y) + threadIdx.y) ;
 453
 454   //zone commune
 455   //forme generique contrib pixel => masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic)
 456   for (ic=0 ; ic<L ; ic++)
 457         {
 458           pix = tex2D(tex_img_inc, j+1, i-1+ic) ;
 459           outval0 += masque[ __umul24(ic,L) +2 ]*pix ;
 460           outval1 += masque[ __umul24(ic,L) +1 ]*pix ;
 461           outval2 += masque[ __umul24(ic,L)    ]*pix ;
 462           pix = tex2D(tex_img_inc, j+2, i-1+ic) ;
 463           outval1 += masque[ __umul24(ic,L) +2 ]*pix ;
 464           outval2 += masque[ __umul24(ic,L) +1 ]*pix ;
 465           outval3 += masque[ __umul24(ic,L)    ]*pix ;
 466           pix = tex2D(tex_img_inc, j+3, i-1+ic) ;
 467           outval2 += masque[ __umul24(ic,L) +2 ]*pix ;
 468           outval3 += masque[ __umul24(ic,L) +1 ]*pix ;
 469           outval4 += masque[ __umul24(ic,L)    ]*pix ;
 470           pix = tex2D(tex_img_inc, j+4, i-1+ic) ;
 471           outval3 += masque[ __umul24(ic,L) +2 ]*pix ;
 472           outval4 += masque[ __umul24(ic,L) +1 ]*pix ;
 473           outval5 += masque[ __umul24(ic,L)    ]*pix ;
 474           pix = tex2D(tex_img_inc, j+5, i-1+ic) ;
 475           outval4 += masque[ __umul24(ic,L) +2 ]*pix ;
 476           outval5 += masque[ __umul24(ic,L) +1 ]*pix ;
 477           outval6 += masque[ __umul24(ic,L)    ]*pix ;
 478           pix = tex2D(tex_img_inc, j+6, i-1+ic) ;
 479           outval5 += masque[ __umul24(ic,L) +2 ]*pix ;
 480           outval6 += masque[ __umul24(ic,L) +1 ]*pix ;
 481           outval7 += masque[ __umul24(ic,L)    ]*pix ;
 482
 483           pix = tex2D(tex_img_inc, j, i-1+ic) ;
 484           outval0 += masque[ __umul24(ic,L) +1 ]*pix ;
 485           outval1 += masque[ __umul24(ic,L)    ]*pix ;
 486           pix = tex2D(tex_img_inc, j-1, i-1+ic) ;
 487           outval0 += masque[ __umul24(ic,L)  ]*pix ;
 488
 489           pix = tex2D(tex_img_inc, j+7, i-1+ic) ;
 490           outval6 += masque[ __umul24(ic,L) +2 ]*pix ;
 491           outval7 += masque[ __umul24(ic,L) +1 ]*pix ;
 492           pix = tex2D(tex_img_inc, j+8, i-1+ic) ;
 493           outval7 += masque[ __umul24(ic,L) +2 ]*pix ;
 494         }
 495
 496   // les 2 pixels hauts
 497   output[ __umul24(i, j_dim) + j   ] = outval0 ;
 498   output[ __umul24(i, j_dim) + j+1 ] = outval1 ;
 499   output[ __umul24(i, j_dim) + j+2 ] = outval2 ;
 500   output[ __umul24(i, j_dim) + j+3 ] = outval3;
 501   output[ __umul24(i, j_dim) + j+4 ] = outval4;
 502   output[ __umul24(i, j_dim) + j+5 ] = outval5;
 503   output[ __umul24(i, j_dim) + j+6 ] = outval6 ;
 504   output[ __umul24(i, j_dim) + j+7 ] = outval7 ;
 505
 506 }
 507
 508
 509
 510
 511
 512 // convolution non separable 8 pixels/thread
 513 // image en texture et noyau en mem constante
 514 __global__ void kernel_convoGene8x8pL5( unsigned char  *output, int j_dim )
 515 {
 516   int ic, base ;
 517   int L=5 ;
 518   unsigned char pix ;
 519   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
 520   float outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
 521
 522   // coordonnees absolues du point de base en haut a gauche
 523   int j = ( __umul24( blockIdx.x, blockDim.x) + threadIdx.x)<< 3 ;
 524   int i = ( __umul24( blockIdx.y, blockDim.y) + threadIdx.y) ;
 525
 526   //zone commune
 527   //forme generique contrib pixel => masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic)
 528   for (ic=0 ; ic<L ; ic++)
 529         {
 530           base = __umul24(ic,L) ;
 531           pix = tex2D(tex_img_inc, j+2, i-2+ic) ;
 532           outval0 += masque[ base +4 ]*pix ;
 533           outval1 += masque[ base +3 ]*pix ;
 534           outval2 += masque[ base +2 ]*pix ;
 535           outval3 += masque[ base +1 ]*pix ;
 536           outval4 += masque[ base    ]*pix ;
 537           pix = tex2D(tex_img_inc, j+3, i-2+ic) ;
 538           outval1 += masque[ base +4 ]*pix ;
 539           outval2 += masque[ base +3 ]*pix ;
 540           outval3 += masque[ base +2 ]*pix ;
 541           outval4 += masque[ base +1 ]*pix ;
 542           outval5 += masque[ base    ]*pix ;
 543           pix = tex2D(tex_img_inc, j+4, i-2+ic) ;
 544           outval2 += masque[ base +4 ]*pix ;
 545           outval3 += masque[ base +3 ]*pix ;
 546           outval4 += masque[ base +2 ]*pix ;
 547           outval5 += masque[ base +1 ]*pix ;
 548           outval6 += masque[ base    ]*pix ;
 549           pix = tex2D(tex_img_inc, j+5, i-2+ic) ;
 550           outval3 += masque[ base +4 ]*pix ;
 551           outval4 += masque[ base +3 ]*pix ;
 552           outval5 += masque[ base +2 ]*pix ;
 553           outval6 += masque[ base +1 ]*pix ;
 554           outval7 += masque[ base    ]*pix ;
 555
 556           pix = tex2D(tex_img_inc, j+1, i-2+ic) ;
 557           outval0 += masque[ base +3 ]*pix ;
 558           outval1 += masque[ base +2 ]*pix ;
 559           outval2 += masque[ base +1 ]*pix ;
 560           outval3 += masque[ base    ]*pix ;
 561
 562           pix = tex2D(tex_img_inc, j, i-2+ic) ;
 563           outval0 += masque[ base +2 ]*pix ;
 564           outval1 += masque[ base +1 ]*pix ;
 565           outval2 += masque[ base    ]*pix ;
 566
 567           pix = tex2D(tex_img_inc, j-1, i-2+ic) ;
 568           outval0 += masque[ base +1 ]*pix ;
 569           outval1 += masque[ base    ]*pix ;
 570
 571           pix = tex2D(tex_img_inc, j-2, i-2+ic) ;
 572           outval0 += masque[ base    ]*pix ;
 573
 574           pix = tex2D(tex_img_inc, j+6, i-2+ic) ;
 575           outval7 += masque[ base +1 ]*pix ;
 576           outval6 += masque[ base +2 ]*pix ;
 577           outval5 += masque[ base +3 ]*pix ;
 578           outval4 += masque[ base +4 ]*pix ;
 579
 580           pix = tex2D(tex_img_inc, j+7, i-2+ic) ;
 581           outval7 += masque[ base +2 ]*pix ;
 582           outval6 += masque[ base +3 ]*pix ;
 583           outval5 += masque[ base +4 ]*pix ;
 584
 585           pix = tex2D(tex_img_inc, j+8, i-2+ic) ;
 586           outval7 += masque[ base +3 ]*pix ;
 587           outval6 += masque[ base +4 ]*pix ;
 588
 589           pix = tex2D(tex_img_inc, j+9, i-2+ic) ;
 590           outval7 += masque[ base +4 ]*pix ;
 591         }
 592
 593   // les 8 pixels
 594   output[ __umul24(i, j_dim) + j   ] = outval0 ;
 595   output[ __umul24(i, j_dim) + j+1 ] = outval1 ;
 596   output[ __umul24(i, j_dim) + j+2 ] = outval2 ;
 597   output[ __umul24(i, j_dim) + j+3 ] = outval3 ;
 598   output[ __umul24(i, j_dim) + j+4 ] = outval4 ;
 599   output[ __umul24(i, j_dim) + j+5 ] = outval5 ;
 600   output[ __umul24(i, j_dim) + j+6 ] = outval6 ;
 601   output[ __umul24(i, j_dim) + j+7 ] = outval7 ;
 602
 603 }
 604
 605
 606 // convolution non separable 4 pixels/thread
 607 // image en texture et noyau en mem constante
 608 __global__ void kernel_convoGene8rx16p( unsigned char  *output, int j_dim, int r )
 609 {
 610   int ic, jc ;
 611   int L=2*r+1, baseIdx ;
 612   unsigned char pix ;
 613   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
 614   float outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
 615   float outval8=0.0, outval9=0.0, outval10=0.0, outval11=0.0 ;
 616   float outval12=0.0, outval13=0.0, outval14=0.0, outval15=0.0 ;
 617
 618   // coordonnees absolues du point de base en haut a gauche
 619   int j = ( __umul24( blockIdx.x, blockDim.x) + threadIdx.x)<< 3 ;
 620   int i = ( __umul24( blockIdx.y, blockDim.y) + threadIdx.y)<< 1 ;
 621
 622   //zone commune
 623   //forme generique contrib pixel => masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic)
 624   for (ic=1 ; ic<L ; ic++)
 625         {
 626           for(jc=1 ; jc<L ; jc++)
 627                 {
 628                   pix = tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
 629                   baseIdx = __umul24(ic-1,L)+jc-1 ;
 630                   outval0 += masque[ baseIdx+L+1 ]*pix ;
 631                   outval1 += masque[ baseIdx+L   ]*pix ;
 632                   outval2 += masque[ baseIdx+1   ]*pix ;
 633                   outval3 += masque[ baseIdx     ]*pix ;
 634
 635                   pix = tex2D(tex_img_inc, j-r+jc+2, i-r+ic) ;
 636                   outval4 += masque[ baseIdx+L+1  ]*pix ;
 637                   outval5 += masque[ baseIdx+L    ]*pix ;
 638                   outval6 += masque[ baseIdx+1    ]*pix ;
 639                   outval7 += masque[ baseIdx      ]*pix ;
 640
 641                   pix = tex2D(tex_img_inc, j-r+jc+4, i-r+ic) ;
 642                   outval8 += masque[ baseIdx+L+1 ]*pix ;
 643                   outval9 += masque[ baseIdx+L   ]*pix ;
 644                   outval10+= masque[ baseIdx+1   ]*pix ;
 645                   outval11+= masque[ baseIdx     ]*pix ;
 646
 647                   pix = tex2D(tex_img_inc, j-r+jc+6, i-r+ic) ;
 648                   outval12 += masque[ baseIdx+L+1 ]*pix ;
 649                   outval13 += masque[ baseIdx+L   ]*pix ;
 650                   outval14 += masque[ baseIdx+1   ]*pix ;
 651                   outval15 += masque[ baseIdx     ]*pix ;
 652
 653                 }
 654         }
 655   // les blocs peripheriques verticaux
 656   for (ic=1 ; ic<L ; ic++)
 657         {
 658           outval0 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r, i-r+ic) ; // bande A pour b0
 659           outval2 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r, i-r+ic) ; // bande A pour b2
 660
 661           outval1 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+1, i-r+ic) ; //bande C pour b1
 662           outval3 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+1, i-r+ic) ; //bande C pour b3
 663
 664           outval4 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r+2, i-r+ic) ; // bande A pour b0
 665           outval6 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r+2, i-r+ic) ; // bande A pour b2
 666
 667           outval5 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+3, i-r+ic) ; //bande C pour b1
 668           outval7 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+3, i-r+ic) ; //bande C pour b3
 669
 670           outval8 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r+4, i-r+ic) ; // bande A pour b0
 671           outval10 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r+4, i-r+ic) ; // bande A pour b2
 672
 673           outval9 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+5, i-r+ic) ; //bande C pour b1
 674           outval11 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+5, i-r+ic) ; //bande C pour b3
 675
 676           outval12 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r+6, i-r+ic) ; // bande A pour b0
 677           outval14 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r+6, i-r+ic) ; // bande A pour b2
 678
 679           outval13 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+7, i-r+ic) ; //bande C pour b1
 680           outval15 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+7, i-r+ic) ; //bande C pour b3
 681         }
 682   // les blocs peripheriques horizontaux
 683   for (jc=1 ; jc<L ; jc++)
 684         {
 685           outval0 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc, i-r) ; // bande B pour b0
 686           outval1 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc, i-r) ; // bande B pour b1
 687
 688           outval2 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc, i+r+1) ; //bande D pour b2
 689           outval3 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc, i+r+1) ; //bande D pour b3
 690
 691           outval4 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc+2, i-r) ; // bande B pour b0
 692           outval5 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc+2, i-r) ; // bande B pour b1
 693
 694           outval6 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc+2, i+r+1) ; //bande D pour b2
 695           outval7 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc+2, i+r+1) ; //bande D pour b3
 696
 697           outval8 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc+4, i-r) ; // bande B pour b0
 698           outval9 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc+4, i-r) ; // bande B pour b1
 699
 700           outval10 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc+4, i+r+1) ; //bande D pour b2
 701           outval11 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc+4, i+r+1) ; //bande D pour b3
 702
 703           outval12 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc+6, i-r) ; // bande B pour b0
 704           outval13 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc+6, i-r) ; // bande B pour b1
 705
 706           outval14 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc+6, i+r+1) ; //bande D pour b2
 707           outval15 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc+6, i+r+1) ; //bande D pour b3
 708         }
 709
 710   //les coins
 711   {
 712         outval0 += masque[ 0                   ]*tex2D(tex_img_inc, j-r  , i-r) ;   // haut gauche pour  b0
 713         outval1 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+1, i-r) ;   // haut droit pour b1
 714         outval2 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r  , i+r+1) ; // bas gauche pour b2
 715         outval3 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+1, i+r+1) ; //bas droit pour b3
 716
 717         outval4 += masque[ 0                   ]*tex2D(tex_img_inc, j-r+2, i-r) ;   // haut gauche pour  b0
 718         outval5 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+3, i-r) ;   // haut droit pour b1
 719         outval6 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r+2, i+r+1) ; // bas gauche pour b2
 720         outval7 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+3, i+r+1) ; //bas droit pour b3
 721
 722         outval8 += masque[ 0                   ]*tex2D(tex_img_inc, j-r+4, i-r) ;   // haut gauche pour  b0
 723         outval9 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+5, i-r) ;   // haut droit pour b1
 724         outval10 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r+4, i+r+1) ; // bas gauche pour b2
 725         outval11 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+5, i+r+1) ; //bas droit pour b3
 726
 727         outval12 += masque[ 0                   ]*tex2D(tex_img_inc, j-r+6, i-r) ;   // haut gauche pour  b0
 728         outval13 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+7, i-r) ;   // haut droit pour b1
 729         outval14 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r+6, i+r+1) ; // bas gauche pour b2
 730         outval15 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+7, i+r+1) ; //bas droit pour b3
 731   }
 732
 733
 734   // les 2 pixels hauts
 735   int base = __umul24(i, j_dim) +j ;
 736   output[ base++   ] = outval0 ;
 737   output[ base++ ] = outval1 ;
 738   output[ base++ ] = outval4 ;
 739   output[ base++ ] = outval5 ;
 740   output[ base++ ] = outval8 ;
 741   output[ base++ ] = outval9 ;
 742   output[ base++ ] = outval12 ;
 743   output[ base++ ] = outval13 ;
 744   // les 2 pixels bas
 745   base = __umul24(i+1, j_dim) +j ;
 746   output[ base++ ] = outval2 ;
 747   output[ base++ ] = outval3 ;
 748   output[ base++ ] = outval6 ;
 749   output[ base++ ] = outval7 ;
 750   output[ base++ ] = outval10 ;
 751   output[ base++ ] = outval11 ;
 752   output[ base++ ] = outval14 ;
 753   output[ base++ ] = outval15 ;
 754
 755 }
 756
 757
 758
 759 // convolution non separable
 760 // image en texture et noyau en mem constante
 761 __global__ void kernel_convoGene8r( unsigned char  *output, int i_dim, int j_dim, int r)
 762 {
 763   int ic, jc ;
 764   int L=2*r+1 ;
 765   float outval0=0.0 ;
 766
 767   // coordonnees absolues du point de base
 768   int j = __mul24( blockIdx.x,blockDim.x ) + threadIdx.x ;
 769   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
 770
 771   for (ic=-r ; ic<=r ; ic++)
 772         for(jc=-r ; jc<=r ; jc++)
 773           outval0 += masque[ __umul24(ic,L)+jc+r ]*tex2D(tex_img_inc, j+jc, i+ic) ;
 774
 775
 776   output[ __umul24(i, j_dim) + j ] = outval0 ;
 777 }
 778
 779 // convolution non separable
 780 // image en texture et noyau en mem constante
 781 __global__ void kernel_convoGeneTex8r( unsigned char  *output, int i_dim, int j_dim, int r)
 782 {
 783   int ic, jc ;
 784   int L=2*r+1 ;
 785   float outval0=0.0 ;
 786
 787   // coordonnees absolues du point de base
 788   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
 789   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
 790
 791   for (ic=0 ; ic<L ; ic++)
 792         for(jc=0 ; jc<L ; jc++)
 793           outval0 += tex1D(tex_noyau, __mul24(ic,L)+jc )*tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
 794
 795   output[ __mul24(i, j_dim) + j ] = outval0 ;
 796 }
 797
 798
 799 __global__ void kernel_convoGene3Reg8( unsigned char  *output, int i_dim, int j_dim )
 800 {
 801   float outval0=0.0 ;
 802   float n0,n1,n2,n3,n4,n5,n6,n7,n8 ;
 803
 804   n0 = (1.0/45) ;
 805   n1 = (2.0/45) ;
 806   n2 = (3.0/45) ;
 807   n3 = (4.0/45) ;
 808   n4 = (5.0/45) ;
 809   n5 = (6.0/45) ;
 810   n6 = (7.0/45) ;
 811   n7 = (8.0/45) ;
 812   n8 = (9.0/45) ;
 813
 814   // coordonnees absolues du point de base
 815   int j = __mul24(blockIdx.x, blockDim.x) + threadIdx.x ;
 816   int i = __mul24(blockIdx.y, blockDim.y) + threadIdx.y ;
 817
 818   outval0 = n8*tex2D(tex_img_inc, j-1, i-1 )
 819               + n7*tex2D(tex_img_inc, j  , i-1 )
 820               + n6*tex2D(tex_img_inc, j+1, i-1 )
 821               + n5*tex2D(tex_img_inc, j-1, i   )
 822               + n4*tex2D(tex_img_inc, j  , i   )
 823           + n3*tex2D(tex_img_inc, j+1, i   )
 824               + n2*tex2D(tex_img_inc, j-1, i+1 )
 825               + n1*tex2D(tex_img_inc, j  , i+1 )
 826               + n0*tex2D(tex_img_inc, j+1, i+1 ) ;
 827
 828   output[ __mul24(i, j_dim) + j  ] = (unsigned char) outval0 ;
 829 }
 830
 831
 832 // convolution non separable
 833 // image en texture et noyau en registres
 834 __global__ void kernel_convoGene5Reg8( unsigned char  *output, int i_dim, int j_dim)
 835 {
 836   float outval0=0.0 ;
 837   float n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n13,n14,n15,n16,n17,n18,n19,n20,n21,n22,n23,n24 ;
 838
 839   n0 = (1.0/25) ;
 840   n1 = (1.0/25) ;
 841   n2 = (1.0/25) ;
 842   n3 = (1.0/25) ;
 843   n4 = (1.0/25) ;
 844   n5 = (1.0/25) ;
 845   n6 = (1.0/25) ;
 846   n7 = (1.0/25) ;
 847   n8 = (1.0/25) ;
 848   n9 = (1.0/25) ;
 849   n10 = (1.0/25) ;
 850   n11 = (1.0/25) ;
 851   n12 = (1.0/25) ;
 852   n13 = (1.0/25) ;
 853   n14 = (1.0/25) ;
 854   n15 = (1.0/25) ;
 855   n16 = (1.0/25) ;
 856   n17 = (1.0/25) ;
 857   n18 = (1.0/25) ;
 858   n19 = (1.0/25) ;
 859   n20 = (1.0/25) ;
 860   n21 = (1.0/25) ;
 861   n22 = (1.0/25) ;
 862   n23 = (1.0/25) ;
 863   n24 = (1.0/25) ;
 864
 865   // coordonnees absolues du point de base
 866   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
 867   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
 868
 869   outval0 = n0*tex2D(tex_img_inc, j-2, i-2 )
 870         + n1*tex2D(tex_img_inc, j-1, i-2 )
 871         + n2*tex2D(tex_img_inc, j  , i-2 )
 872         + n3*tex2D(tex_img_inc, j+1, i-2 )
 873         + n4*tex2D(tex_img_inc, j+2, i-2 )
 874     + n5*tex2D(tex_img_inc, j-2, i-1 )
 875         + n6*tex2D(tex_img_inc, j-1, i-1 )
 876         + n7*tex2D(tex_img_inc, j  , i-1 )
 877         + n8*tex2D(tex_img_inc, j+1, i-1 )
 878         + n9*tex2D(tex_img_inc, j+2, i-1 )
 879   + n10*tex2D(tex_img_inc, j-2, i )
 880   + n11*tex2D(tex_img_inc, j-1, i )
 881   + n12*tex2D(tex_img_inc, j  , i )
 882   + n13*tex2D(tex_img_inc, j+1, i )
 883   + n14*tex2D(tex_img_inc, j+2, i )
 884         + n15*tex2D(tex_img_inc, j-2, i+1 )
 885         + n16*tex2D(tex_img_inc, j-1, i+1 )
 886         + n17*tex2D(tex_img_inc, j  , i+1 )
 887         + n18*tex2D(tex_img_inc, j+1, i+1 )
 888         + n19*tex2D(tex_img_inc, j+2, i+1 )
 889         + n20*tex2D(tex_img_inc, j-2, i+2 )
 890   + n21*tex2D(tex_img_inc, j-1, i+2 )
 891   + n22*tex2D(tex_img_inc, j  , i+2 )
 892   + n23*tex2D(tex_img_inc, j+1, i+2 )
 893         + n24*tex2D(tex_img_inc, j+2, i+2 ) ;
 894
 895   output[ __mul24(i, j_dim) + j  ] = outval0 ;
 896 }
 897
 898 // convolution non separable
 899 // image en texture et noyau en registres
 900 __global__ void kernel_convoGene7Reg8( unsigned char  *output, int i_dim, int j_dim)
 901 {
 902   float outval0=0.0 ;
 903   float n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n13,n14,n15,n16,n17,n18,n19,n20,n21,n22,n23,n24 ;
 904   float n25,n26,n27,n28,n29,n30,n31,n32,n33,n34,n35,n36,n37,n38,n39,n40,n41,n42,n43,n44,n45,n46,n47,n48 ;
 905
 906   n0 = (1.0/49) ;
 907   n1 = (1.0/49) ;
 908   n2 = (1.0/49) ;
 909   n3 = (1.0/49) ;
 910   n4 = (1.0/49) ;
 911   n5 = (1.0/49) ;
 912   n6 = (1.0/49) ;
 913   n7 = (1.0/49) ;
 914   n8 = (1.0/49) ;
 915   n9 = (1.0/49) ;
 916   n10 = (1.0/49) ;
 917   n11 = (1.0/49) ;
 918   n12 = (1.0/49) ;
 919   n13 = (1.0/49) ;
 920   n14 = (1.0/49) ;
 921   n15 = (1.0/49) ;
 922   n16 = (1.0/49) ;
 923   n17 = (1.0/49) ;
 924   n18 = (1.0/49) ;
 925   n19 = (1.0/49) ;
 926   n20 = (1.0/49) ;
 927   n21 = (1.0/49) ;
 928   n22 = (1.0/49) ;
 929   n23 = (1.0/49) ;
 930   n24 = (1.0/49) ;
 931
 932   n25 = (1.0/49) ;
 933   n26 = (1.0/49) ;
 934   n27 = (1.0/49) ;
 935   n28 = (1.0/49) ;
 936   n29 = (1.0/49) ;
 937   n30 = (1.0/49) ;
 938   n31 = (1.0/49) ;
 939   n32 = (1.0/49) ;
 940   n33 = (1.0/49) ;
 941   n34 = (1.0/49) ;
 942   n35 = (1.0/49) ;
 943   n36 = (1.0/49) ;
 944   n37 = (1.0/49) ;
 945   n38 = (1.0/49) ;
 946   n39 = (1.0/49) ;
 947   n40 = (1.0/49) ;
 948   n41 = (1.0/49) ;
 949   n42 = (1.0/49) ;
 950   n43 = (1.0/49) ;
 951   n44 = (1.0/49) ;
 952   n45 = (1.0/49) ;
 953   n46 = (1.0/49) ;
 954   n47 = (1.0/49) ;
 955   n48 = (1.0/49) ;
 956
 957
 958
 959   // coordonnees absolues du point de base
 960   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
 961   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
 962
 963   outval0 =
 964           n0*tex2D(tex_img_inc, j-3, i-3 )
 965         + n1*tex2D(tex_img_inc, j-2, i-3 )
 966         + n2*tex2D(tex_img_inc, j-1, i-3 )
 967         + n3*tex2D(tex_img_inc, j  , i-3 )
 968         + n4*tex2D(tex_img_inc, j+1, i-3 )
 969     + n5*tex2D(tex_img_inc, j+2, i-3 )
 970         + n6*tex2D(tex_img_inc, j+3, i-3 )
 971         +  n7*tex2D(tex_img_inc, j-3, i-2 )
 972         +  n8*tex2D(tex_img_inc, j-2, i-2 )
 973         +  n9*tex2D(tex_img_inc, j-1, i-2 )
 974         + n10*tex2D(tex_img_inc, j  , i-2 )
 975         + n11*tex2D(tex_img_inc, j+1, i-2 )
 976     + n12*tex2D(tex_img_inc, j+2, i-2 )
 977         + n13*tex2D(tex_img_inc, j+3, i-2 )
 978         + n14*tex2D(tex_img_inc, j-3, i-1 )
 979         + n15*tex2D(tex_img_inc, j-2, i-1 )
 980         + n16*tex2D(tex_img_inc, j-1, i-1 )
 981         + n17*tex2D(tex_img_inc, j  , i-1 )
 982         + n18*tex2D(tex_img_inc, j+1, i-1 )
 983     + n19*tex2D(tex_img_inc, j+2, i-1 )
 984         + n20*tex2D(tex_img_inc, j+3, i-1 )
 985         + n21*tex2D(tex_img_inc, j-3, i )
 986         + n22*tex2D(tex_img_inc, j-2, i )
 987         + n23*tex2D(tex_img_inc, j-1, i )
 988         + n24*tex2D(tex_img_inc, j  , i )
 989         + n25*tex2D(tex_img_inc, j+1, i )
 990     + n26*tex2D(tex_img_inc, j+2, i )
 991         + n27*tex2D(tex_img_inc, j+3, i )
 992         + n28*tex2D(tex_img_inc, j-3, i+1 )
 993         + n29*tex2D(tex_img_inc, j-2, i+1 )
 994         + n30*tex2D(tex_img_inc, j-1, i+1 )
 995         + n31*tex2D(tex_img_inc, j  , i+1 )
 996         + n32*tex2D(tex_img_inc, j+1, i+1 )
 997     + n33*tex2D(tex_img_inc, j+2, i+1 )
 998         + n34*tex2D(tex_img_inc, j+3, i+1 )
 999         + n35*tex2D(tex_img_inc, j-3, i+2 )
1000         + n36*tex2D(tex_img_inc, j-2, i+2 )
1001         + n37*tex2D(tex_img_inc, j-1, i+2 )
1002         + n38*tex2D(tex_img_inc, j  , i+2 )
1003         + n39*tex2D(tex_img_inc, j+1, i+2 )
1004     + n40*tex2D(tex_img_inc, j+2, i+2 )
1005         + n41*tex2D(tex_img_inc, j+3, i+2 )
1006         + n42*tex2D(tex_img_inc, j-3, i+3 )
1007         + n43*tex2D(tex_img_inc, j-2, i+3 )
1008         + n44*tex2D(tex_img_inc, j-1, i+3 )
1009         + n45*tex2D(tex_img_inc, j  , i+3 )
1010         + n46*tex2D(tex_img_inc, j+1, i+3 )
1011     + n47*tex2D(tex_img_inc, j+2, i+3 )
1012         + n48*tex2D(tex_img_inc, j+3, i+3 ) ;
1013
1014
1015   output[ __mul24(i, j_dim) + j  ] = outval0 ;
1016 }
1017
1018
1019 __global__ void kernel_convoNonSep16( unsigned short*output, int i_dim, int j_dim)
1020 {
1021   int idb, ic, jc ;
1022   int r = rnoyau ;
1023   int L = 2*r+1 ;
1024   int N = L*L ;
1025   float outval0=0 ;
1026
1027   // coordonnees absolues du point de base
1028   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
1029   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1030
1031   for (idb=0 ; idb< N ; idb++)
1032         {
1033           ic = i-r + idb/L ;
1034           jc = j-r + idb - (idb/L)*L ;
1035           outval0 += ( noyau[ idb ]*tex2D(tex_img_in, jc, ic)) ;
1036         }
1037
1038   output[ __mul24(i, j_dim) + j  ] = outval0 ;
1039 }
1040
1041 // convolution non separable
1042 // image en texture et noyau en mem constante
1043 // prefetch des pixels en smem
1044 // 1 pixel par thread
1045 __global__ void kernel_convoNonSepSh( unsigned short*output, int i_dim, int j_dim)
1046 {
1047   int idb, ic, jc ;
1048   int L = 2*rnoyau+1 ;
1049   float outval0=0 ;
1050
1051   // coordonnees absolues du point de base
1052   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
1053   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1054   int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
1055
1056   // chargement en smem
1057   int idrow = threadIdx.y*(blockDim.x+2*rnoyau) ;
1058
1059   extern __shared__ int roi[];
1060
1061   // bloc 0 (en haut à gauche)
1062   roi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-rnoyau, i-rnoyau) ;
1063   // bloc 1 (en haut à droite)...
1064   if ( threadIdx.x < L-1 ) //...ou plutot ce qu'il en manque
1065         roi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-rnoyau, i-rnoyau ) ;
1066   // bloc 2 ( en bas à gauche)
1067   if ( threadIdx.y < L-1 )
1068         {
1069           idrow = (threadIdx.y+blockDim.y)*(blockDim.x+L-1) ;
1070           roi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-rnoyau, i+blockDim.y-rnoyau ) ;
1071           //bloc 4 ( en bas à doite )...
1072           if ( threadIdx.x < L-1 ) //...ou ce qu'il en manque
1073                 roi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-rnoyau, i+blockDim.y-rnoyau ) ;
1074         }
1075   __syncthreads();
1076
1077   // calculs de convolution
1078   for (ic=0 ; ic<L ; ic++)
1079         for( jc=0 ; jc<L ; jc++)
1080           outval0 += ( noyau[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(blockDim.x+L-1)) + jc+threadIdx.x]) ;
1081
1082   // 1 pixel par thread --> global mem
1083   output[ idx  ] = outval0 ;
1084
1085 }
1086
1087 __global__ void kernel_convoNonSepSh_2p(unsigned short *output, int i_dim, int j_dim)
1088 {
1089   int idb, ic, jc ;
1090   int L = 2*rnoyau+1 ;
1091   float outval0=0.0, outval1=0.0 ;
1092   int bdimX = blockDim.x<<1 ;
1093   int tidX = threadIdx.x<<1 ;
1094
1095   // coordonnees absolues du point de base
1096   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
1097   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1098   int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
1099
1100   // chargement en smem
1101   int idrow = threadIdx.y*(bdimX+L-1) ;
1102
1103   extern __shared__ int roi[];
1104
1105   // bloc 0 (en haut à gauche)
1106   roi[ idrow  + tidX ] = tex2D(tex_img_ins, j-rnoyau, i-rnoyau) ;
1107   roi[ idrow  + tidX +1] = tex2D(tex_img_ins, j-rnoyau+1, i-rnoyau) ;
1108
1109   // bloc 1 (en haut à droite)...
1110   if ( threadIdx.x < rnoyau ) //...ou plutot ce qu'il en manque
1111         {
1112           roi[ idrow + tidX+bdimX   ] = tex2D( tex_img_ins, j+bdimX-rnoyau, i-rnoyau ) ;
1113           roi[ idrow + tidX+bdimX +1] = tex2D( tex_img_ins, j+bdimX-rnoyau+1, i-rnoyau ) ;
1114         }
1115
1116   // bloc 2 ( en bas à gauche)
1117   if ( threadIdx.y < L-1 )
1118         {
1119           idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
1120           roi[ idrow + tidX ] = tex2D( tex_img_ins, j-rnoyau, i+blockDim.y-rnoyau ) ;
1121           roi[ idrow + tidX +1] = tex2D( tex_img_ins, j-rnoyau+1, i+blockDim.y-rnoyau ) ;
1122           //bloc 4 ( en bas à doite )...
1123
1124           if ( 2*threadIdx.x < L-1 ) //...ou ce qu'il en manque
1125                 {
1126                   roi[ idrow + tidX+bdimX ] = tex2D( tex_img_ins, j+bdimX-rnoyau, i+blockDim.y-rnoyau ) ;
1127                   roi[ idrow + tidX+bdimX +1] = tex2D( tex_img_ins, j+bdimX-rnoyau+1, i+blockDim.y-rnoyau ) ;
1128                   }
1129         }
1130   __syncthreads();
1131
1132   // calculs de convolution
1133   for (ic=0 ; ic<L ; ic++)
1134         for( jc=0 ; jc<L ; jc++)
1135           {
1136                 outval0 += (noyau[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(2*blockDim.x+L-1)) + jc+2*threadIdx.x]) ;
1137                 outval1 += (noyau[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(2*blockDim.x+L-1)) + jc+2*threadIdx.x+1]) ;
1138           }
1139
1140   // 1 pixel par thread --> global mem
1141   output[ idx  ] = outval0 ;
1142   output[ idx+1 ] = outval1 ;
1143 }
1144
1145 __global__ void kernel_convoNonSepSh_2p(unsigned char *output, int i_dim, int j_dim, int r)
1146 {
1147   int idb, ic, jc ;
1148   int L = 2*r+1 ;
1149   float outval0=0.0, outval1=0.0 ;
1150   int bdimX = blockDim.x<<1 ;
1151   int tidX = threadIdx.x<<1 ;
1152
1153   // coordonnees absolues du point de base
1154   int j = __umul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
1155   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1156   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
1157
1158   // chargement en smem
1159   int idrow = threadIdx.y*(bdimX+L-1) ;
1160
1161   extern __shared__ int roi[];
1162
1163   // bloc 0 (en haut à gauche)
1164   roi[ idrow  + tidX    ] = tex2D( tex_img_inc, j-r  , i-r) ;
1165   roi[ idrow  + tidX +1 ] = tex2D( tex_img_inc, j-r+1, i-r) ;
1166
1167   // bloc 1 (en haut à droite)...
1168   if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
1169         {
1170           roi[ idrow + bdimX + tidX    ] = tex2D( tex_img_inc, j+bdimX-r  , i-r ) ;
1171           roi[ idrow + bdimX + tidX +1 ] = tex2D( tex_img_inc, j+bdimX-r+1, i-r ) ;
1172           }
1173
1174   // bloc 2 ( en bas à gauche)
1175   if ( threadIdx.y < L-1 )
1176         {
1177           idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
1178           roi[ idrow + tidX ] = tex2D( tex_img_inc, j-r, i+blockDim.y-r ) ;
1179           roi[ idrow + tidX +1] = tex2D( tex_img_inc, j-r+1, i+blockDim.y-r ) ;
1180           //bloc 4 ( en bas à doite )...
1181
1182           if ( threadIdx.x < r ) //...ou ce qu'il en manque
1183                 {
1184                   roi[ idrow + tidX+bdimX   ] = tex2D( tex_img_inc, j+bdimX-r  , i+blockDim.y-r ) ;
1185                   roi[ idrow + tidX+bdimX +1] = tex2D( tex_img_inc, j+bdimX-r+1, i+blockDim.y-r ) ;
1186                   }
1187         }
1188   __syncthreads();
1189
1190   // calculs de convolution
1191   for (ic=0 ; ic<L ; ic++)
1192         for( jc=0 ; jc<L ; jc++)
1193           {
1194                 outval0 += (masque[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX]) ;
1195                 outval1 += (masque[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX+1]) ;
1196           }
1197
1198   // 1 pixel par thread --> global mem
1199   output[ idx  ] = outval0 ;
1200   output[ idx+1 ] = outval1 ;
1201 }
1202
1203 //4 pixels en ligne par thread
1204 __global__ void kernel_convoNonSepSh_4p(unsigned char *output, int i_dim, int j_dim, int r)
1205 {
1206   int idb, ic, jc ;
1207   int L = 2*r+1 ;
1208   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
1209   int bdimX = blockDim.x<<2 ;
1210   int tidX = threadIdx.x<<2 ;
1211
1212   // coordonnees absolues du point de base
1213   int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<2 ;
1214   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1215   int j0= __umul24(blockIdx.x,blockDim.x)<<2 ;
1216
1217   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
1218
1219
1220   // chargement en smem
1221   int idrow = threadIdx.y*(bdimX+L-1) ;
1222
1223   extern __shared__ unsigned char roi4p[];
1224
1225   // bloc 0 (en haut à gauche)
1226   roi4p[ idrow  + tidX   ] = tex2D(tex_img_inc, j-r  , i-r) ;
1227   roi4p[ idrow  + tidX +1] = tex2D(tex_img_inc, j-r+1, i-r) ;
1228   roi4p[ idrow  + tidX +2] = tex2D(tex_img_inc, j-r+2, i-r) ;
1229   roi4p[ idrow  + tidX +3] = tex2D(tex_img_inc, j-r+3, i-r) ;
1230
1231   // bloc 1 (en haut à droite)...
1232   if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
1233         {
1234           roi4p[ idrow + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r ) ;
1235           roi4p[ idrow + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r ) ;
1236         }
1237   // bloc 2 ( en bas à gauche)
1238   if ( threadIdx.y < L-1 )
1239         {
1240           idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
1241           roi4p[ idrow + tidX   ] = tex2D( tex_img_inc, j-r  , i+blockDim.y-r ) ;
1242           roi4p[ idrow + tidX +1] = tex2D( tex_img_inc, j-r+1, i+blockDim.y-r ) ;
1243           roi4p[ idrow + tidX +2] = tex2D( tex_img_inc, j-r+2, i+blockDim.y-r ) ;
1244           roi4p[ idrow + tidX +3] = tex2D( tex_img_inc, j-r+3, i+blockDim.y-r ) ;
1245
1246           //bloc 4 ( en bas à doite )...
1247           if ( threadIdx.x < r ) //...ou ce qu'il en manque
1248                 {
1249                   roi4p[ idrow + bdimX +threadIdx.x   ] = tex2D( tex_img_inc, j0-r +bdimX +threadIdx.x, i+blockDim.y-r ) ;
1250                   roi4p[ idrow + bdimX +threadIdx.x +r] = tex2D( tex_img_inc, j0   +bdimX +threadIdx.x, i+blockDim.y-r ) ;
1251                 }
1252         }
1253   __syncthreads();
1254
1255   // calculs de convolution
1256   for (ic=0 ; ic<L ; ic++)
1257         for( jc=0 ; jc<L ; jc++)
1258           {
1259                 int baseRoi = __umul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX ;
1260                 float valMask = masque[ __umul24(ic,L)+jc ] ;
1261
1262                 outval0 += valMask*roi4p[ baseRoi ] ;
1263                 outval1 += valMask*roi4p[ baseRoi +1 ] ;
1264                 outval2 += valMask*roi4p[ baseRoi +2 ] ;
1265                 outval3 += valMask*roi4p[ baseRoi +3 ] ;
1266           }
1267
1268   // 1 pixel par thread --> global mem
1269   output[ idx   ] = outval0 ;
1270   output[ idx+1 ] = outval1 ;
1271   output[ idx+2 ] = outval2 ;
1272   output[ idx+3 ] = outval3 ;
1273 }
1274
1275
1276 //4 pixels en carre par thread
1277 __global__ void kernel_convoNonSepSh_4pcarre(unsigned char *output, int i_dim, int j_dim, int r)
1278 {
1279   int idb, ic, jc ;
1280   int L = 2*r+1 ;
1281   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
1282   int bdimX = blockDim.x<<1 ;
1283   int tidX = threadIdx.x<<1 ;
1284   int bdimY = blockDim.y<<1 ;
1285   int tidY = threadIdx.y<<1 ;
1286
1287   // coordonnees absolues du point de base
1288   int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<1 ;
1289   int i = (__umul24( blockIdx.y, blockDim.y) + threadIdx.y)<<1 ;
1290   int j0= __umul24(blockIdx.x,blockDim.x)<<1 ;
1291   int i0= __umul24(blockIdx.y,blockDim.y)<<1 ;
1292
1293   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
1294
1295
1296   // chargement en smem
1297   int idrowBase   = tidY*(bdimX+L-1) ;
1298   int idrowSecond = (tidY+1)*(bdimX+L-1) ;
1299
1300
1301   extern __shared__ unsigned char roi4p2[];
1302
1303   // bloc 0 (en haut à gauche)
1304   roi4p2[ idrowBase  + tidX   ] = tex2D(tex_img_inc, j-r  , i-r  ) ;
1305   roi4p2[ idrowBase  + tidX +1] = tex2D(tex_img_inc, j-r+1, i-r  ) ;
1306   roi4p2[ idrowSecond+ tidX   ] = tex2D(tex_img_inc, j-r  , i-r+1) ;
1307   roi4p2[ idrowSecond+ tidX +1] = tex2D(tex_img_inc, j-r+1, i-r+1) ;
1308
1309   // bloc 1 (en haut à droite)...
1310   if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
1311         {
1312           roi4p2[ idrowBase + bdimX   + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r ) ;
1313           roi4p2[ idrowBase + bdimX   + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r ) ;
1314           roi4p2[ idrowSecond + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r+1 ) ;
1315           roi4p2[ idrowSecond + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r+1 ) ;
1316         }
1317   // bloc 2 ( en bas à gauche)
1318   if ( threadIdx.y < L-1 )
1319         {
1320           idrowBase   = (bdimY + threadIdx.y)*(bdimX+L-1) ;
1321           roi4p2[ idrowBase + tidX     ] = tex2D( tex_img_inc, j-r  , i0-r +bdimY +threadIdx.y ) ;
1322           roi4p2[ idrowBase + tidX   +1] = tex2D( tex_img_inc, j-r+1, i0-r +bdimY +threadIdx.y ) ;
1323
1324           //bloc 4 ( en bas à doite )...
1325           if ( threadIdx.x < r ) //...ou ce qu'il en manque
1326                 {
1327                   roi4p2[ idrowBase + bdimX + threadIdx.x    ]   = tex2D( tex_img_inc, j0-r +bdimX +threadIdx.x  , i0-r +bdimY +threadIdx.y ) ;
1328                   roi4p2[ idrowBase + bdimX + threadIdx.x +r ]   = tex2D( tex_img_inc, j0   +bdimX +threadIdx.x  , i0-r +bdimY +threadIdx.y ) ;
1329                 }
1330         }
1331   __syncthreads();
1332
1333   // calculs de convolution
1334   for (ic=0 ; ic<L ; ic++)
1335         for( jc=0 ; jc<L ; jc++)
1336           {
1337                 int baseRoi = __umul24(ic+tidY,(bdimX+L-1)) + jc+tidX ;
1338                 float valMask = masque[ __umul24(ic,L)+jc ] ;
1339
1340                 outval0 += valMask*roi4p2[ baseRoi ] ;
1341                 outval1 += valMask*roi4p2[ baseRoi +1 ] ;
1342                 outval2 += valMask*roi4p2[ baseRoi +bdimX+L-1 ] ;
1343                 outval3 += valMask*roi4p2[ baseRoi +bdimX+L ] ;
1344           }
1345
1346   // 1 pixel par thread --> global mem
1347   output[ idx   ] = outval0 ;
1348   output[ idx+1 ] = outval1 ;
1349   output[ idx+j_dim ] = outval2 ;
1350   output[ idx+j_dim+1 ] = outval3 ;
1351 }
1352
1353
1354 //8 pixels en ligne par thread
1355 __global__ void kernel_convoNonSepSh_8p(unsigned char *output, int i_dim, int j_dim, int r)
1356 {
1357   int idb, ic, jc;
1358   int L = 2*r+1 ;
1359   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
1360   int bdimX = blockDim.x<<3 ;
1361   int tidX = threadIdx.x<<3 ;
1362
1363   // coordonnees absolues du point de base
1364   int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<3 ;
1365   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1366   int j0= __umul24(blockIdx.x,blockDim.x)<<3 ;
1367
1368   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
1369
1370
1371   // chargement en smem
1372   int idrow = threadIdx.y*(bdimX+L-1) ;
1373
1374   extern __shared__ unsigned char roi8p[];
1375
1376   // bloc 0 (en haut à gauche)
1377   for (int p=0; p<8; p++)
1378   roi8p[ idrow  + tidX +p ] = tex2D(tex_img_inc, j-r+p  , i-r) ;
1379
1380   // bloc 1 (en haut à droite)...
1381   if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
1382         {
1383           roi8p[ idrow + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r ) ;
1384           roi8p[ idrow + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r ) ;
1385         }
1386   // bloc 2 ( en bas à gauche)
1387   if ( threadIdx.y < L-1 )
1388         {
1389           idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
1390           for (int p=0; p<8; p++)
1391                 roi8p[ idrow + tidX +p  ] = tex2D( tex_img_inc, j-r+p  , i+blockDim.y-r ) ;
1392
1393           //bloc 4 ( en bas à doite )...
1394           if ( threadIdx.x < r ) //...ou ce qu'il en manque
1395                 {
1396                   roi8p[ idrow + bdimX +threadIdx.x   ] = tex2D( tex_img_inc, j0-r +bdimX +threadIdx.x, i+blockDim.y-r ) ;
1397                   roi8p[ idrow + bdimX +threadIdx.x +r] = tex2D( tex_img_inc, j0   +bdimX +threadIdx.x, i+blockDim.y-r ) ;
1398                 }
1399         }
1400   __syncthreads();
1401
1402   // calculs de convolution
1403   for (ic=0 ; ic<L ; ic++)
1404         for( jc=0 ; jc<L ; jc++)
1405           {
1406                 int baseRoi = __umul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX ;
1407                 float valMask = masque[ __umul24(ic,L)+jc ] ;
1408                 outval0 += valMask*roi8p[ baseRoi ] ;
1409                 outval1 += valMask*roi8p[ baseRoi +1 ] ;
1410                 outval2 += valMask*roi8p[ baseRoi +2 ] ;
1411                 outval3 += valMask*roi8p[ baseRoi +3 ] ;
1412                 outval4 += valMask*roi8p[ baseRoi +4 ] ;
1413                 outval5 += valMask*roi8p[ baseRoi +5 ] ;
1414                 outval6 += valMask*roi8p[ baseRoi +6 ] ;
1415                 outval7 += valMask*roi8p[ baseRoi +7 ] ;
1416           }
1417
1418   // 1 pixel par thread --> global mem
1419   output[ idx   ] = outval0 ;
1420   output[ idx+1 ] = outval1 ;
1421   output[ idx+2 ] = outval2 ;
1422   output[ idx+3 ] = outval3 ;
1423   output[ idx+4 ] = outval4 ;
1424   output[ idx+5 ] = outval5 ;
1425   output[ idx+6 ] = outval6 ;
1426   output[ idx+7 ] = outval7 ;
1427 }
1428
1429
1430 //16 pixels en ligne par thread
1431 __global__ void kernel_convoNonSepSh_16p(unsigned char *output, int i_dim, int j_dim, int r)
1432 {
1433   int idb, ic, jc;
1434   int L = 2*r+1 ;
1435   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
1436   float outval8=0.0, outval9=0.0, outval10=0.0, outval11=0.0, outval12=0.0, outval13=0.0, outval14=0.0, outval15=0.0 ;
1437   int bdimX = blockDim.x<<4 ;
1438   int tidX = threadIdx.x<<4 ;
1439
1440   // coordonnees absolues du point de base
1441   int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<4 ;
1442   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1443   int j0= __umul24(blockIdx.x,blockDim.x)<<4 ;
1444
1445   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
1446
1447
1448   // chargement en smem
1449   int idrow = threadIdx.y*(bdimX+L-1) ;
1450
1451   extern __shared__ unsigned char roi16p[];
1452
1453   // bloc 0 (en haut à gauche)
1454   for (int p=0; p<16; p++)
1455         roi16p[ idrow  + tidX +p ] = tex2D(tex_img_inc, j-r+p  , i-r) ;
1456
1457   // bloc 1 (en haut à droite)...
1458   if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
1459         {
1460           roi16p[ idrow + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r ) ;
1461           roi16p[ idrow + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r ) ;
1462         }
1463   // bloc 2 ( en bas à gauche)
1464   if ( threadIdx.y < L-1 )
1465         {
1466           idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
1467           for (int p=0; p<16; p++)
1468                 roi16p[ idrow + tidX +p  ] = tex2D( tex_img_inc, j-r+p  , i+blockDim.y-r ) ;
1469
1470           //bloc 4 ( en bas à doite )...
1471           if ( threadIdx.x < r ) //...ou ce qu'il en manque
1472                 {
1473                   roi16p[ idrow + bdimX +threadIdx.x   ] = tex2D( tex_img_inc, j0-r +bdimX +threadIdx.x, i+blockDim.y-r ) ;
1474                   roi16p[ idrow + bdimX +threadIdx.x +r] = tex2D( tex_img_inc, j0   +bdimX +threadIdx.x, i+blockDim.y-r ) ;
1475                 }
1476         }
1477   __syncthreads();
1478
1479   // calculs de convolution
1480   for (ic=0 ; ic<L ; ic++)
1481         for( jc=0 ; jc<L ; jc++)
1482           {
1483                 int baseRoi = __umul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX ;
1484                 float valMask = masque[ __umul24(ic,L)+jc ] ;
1485                 outval0 += valMask*roi16p[ baseRoi ] ;
1486                 outval1 += valMask*roi16p[ baseRoi +1 ] ;
1487                 outval2 += valMask*roi16p[ baseRoi +2 ] ;
1488                 outval3 += valMask*roi16p[ baseRoi +3 ] ;
1489                 outval4 += valMask*roi16p[ baseRoi +4 ] ;
1490                 outval5 += valMask*roi16p[ baseRoi +5 ] ;
1491                 outval6 += valMask*roi16p[ baseRoi +6 ] ;
1492                 outval7 += valMask*roi16p[ baseRoi +7 ] ;
1493                 outval8 += valMask*roi16p[ baseRoi +8] ;
1494                 outval9 += valMask*roi16p[ baseRoi +9 ] ;
1495                 outval10 += valMask*roi16p[ baseRoi +10 ] ;
1496                 outval11 += valMask*roi16p[ baseRoi +11 ] ;
1497                 outval12 += valMask*roi16p[ baseRoi +12 ] ;
1498                 outval13 += valMask*roi16p[ baseRoi +13 ] ;
1499                 outval14 += valMask*roi16p[ baseRoi +14 ] ;
1500                 outval15 += valMask*roi16p[ baseRoi +15 ] ;
1501           }
1502
1503   // 1 pixel par thread --> global mem
1504   output[ idx   ] = outval0 ;
1505   output[ idx+1 ] = outval1 ;
1506   output[ idx+2 ] = outval2 ;
1507   output[ idx+3 ] = outval3 ;
1508   output[ idx+4 ] = outval4 ;
1509   output[ idx+5 ] = outval5 ;
1510   output[ idx+6 ] = outval6 ;
1511   output[ idx+7 ] = outval7 ;
1512   output[ idx+8 ] = outval8 ;
1513   output[ idx+9 ] = outval9 ;
1514   output[ idx+10 ] = outval10 ;
1515   output[ idx+11 ] = outval11 ;
1516   output[ idx+12 ] = outval12 ;
1517   output[ idx+13 ] = outval13 ;
1518   output[ idx+14 ] = outval14 ;
1519   output[ idx+15 ] = outval15 ;
1520 }
1521
1522
1523
1524 // convolution non separable
1525 // image en texture et noyau en mem constante
1526 // fetch direct des pixels
1527 // 2 pixels traités par thread => meilleur débit
1528 __global__ void kernel_convoNonSep_2p( int *output, int i_dim, int j_dim)
1529 {
1530   int idb, ic, jc ;
1531   int r = rnoyau ;
1532   int L = 2*r+1 ;
1533   int N = L*L ;
1534   float outval0=0, outval1=0 ;
1535
1536   // coordonnees absolues du point de base
1537   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x, 2) ;
1538   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1539   int idx = __mul24(i, j_dim) + j ;              // indice  dans l'image
1540
1541   #pragma unroll
1542   for (idb=0 ; idb< N ; idb++)
1543         {
1544           ic = i-r + idb/L ;
1545           jc = j-r + idb - (idb/L)*L ;
1546           outval0 += ( noyau[ idb ]*tex2D(tex_img_in, jc, ic)) ;
1547           outval1 += ( noyau[ idb ]*tex2D(tex_img_in, jc+1, ic)) ;
1548         }
1549
1550   output[ idx   ] = outval0  ;
1551   output[ idx+1 ] = outval1  ;
1552
1553
1554 }
1555
1556
1557 // convolution non separable
1558 // image en texture et noyau en mem constante
1559 // fetch direct des pixels
1560 // 2 pixels traités par thread => meilleur débit
1561 __global__ void kernel_convoNonSep_2p_s( unsigned short*output, int i_dim, int j_dim)
1562 {
1563   int idb, ic, jc ;
1564   int r = rnoyau ;
1565   int L = 2*r+1 ;
1566   int N = L*L ;
1567   float outval0=0, outval1=0 ;
1568
1569   // coordonnees absolues du point de base
1570   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x, 2) ;
1571   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1572   int idx = __mul24(i, j_dim) + j ;              // indice  dans l'image
1573
1574   #pragma unroll
1575   for (idb=0 ; idb< N ; idb++)
1576         {
1577           ic = i-r + idb/L ;
1578           jc = j-r + idb - (idb/L)*L ;
1579           outval0 += ( noyau[ idb ]*tex2D(tex_img_ins, jc, ic)) ;
1580           outval1 += ( noyau[ idb ]*tex2D(tex_img_ins, jc+1, ic)) ;
1581         }
1582
1583   output[ idx   ] = outval0  ;
1584   output[ idx+1 ] = outval1  ;
1585
1586 }
1587
1588
1589
1590 // convolution separable
1591 // image en texture et noyau en mem constante
1592 __global__ void kernel_convoSep8V( unsigned char  *output, int i_dim, int j_dim, int r)
1593 {
1594   int ic ;
1595   int L=2*r+1 ;
1596   float outval0=0.0 ;
1597
1598   // coordonnees absolues du point de base
1599   int j = __mul24( blockIdx.x, blockDim.x ) + threadIdx.x ;
1600   int i = __mul24( blockIdx.y, blockDim.y ) + threadIdx.y ;
1601
1602   //vertical
1603   for (ic=0 ; ic<L ; ic++)
1604           outval0 += masque[ ic ]*tex2D(tex_img_inc, j, i-r+ic) ;
1605
1606   output[ __mul24(i, j_dim) + j ] = outval0 ;
1607 }
1608
1609
1610 __global__ void kernel_convoSep8H( unsigned char  *output, int i_dim, int j_dim, int r)
1611 {
1612   int jc ;
1613   int L=2*r+1 ;
1614   float outval0=0.0 ;
1615
1616   // coordonnees absolues du point de base
1617   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
1618   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1619
1620   //horizontal
1621   for (jc=0 ; jc<L ; jc++)
1622           outval0 += masque[ L+jc ]*tex2D(tex_img_inc, j-r+jc, i) ;
1623
1624   output[ __mul24(i, j_dim) + j ] = outval0 ;
1625 }
1626
1627 // convolution separable
1628 // image en texture et noyau en mem constante
1629 __global__ void kernel_convoSep8Vx8p( unsigned char  *output, int i_dim, int j_dim, int r)
1630 {
1631   int ic, baseIdx ;
1632   int L=2*r+1 ;
1633   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
1634
1635   // coordonnees absolues du point de base
1636   int j = (__mul24( blockIdx.x, blockDim.x ) + threadIdx.x)<<3 ;
1637   int i = __mul24( blockIdx.y, blockDim.y ) + threadIdx.y ;
1638
1639   //vertical
1640   for (ic=0 ; ic<L ; ic++)
1641         {
1642           float valMask = masque[ ic ] ;
1643           outval0 += valMask*tex2D(tex_img_inc, j  , i-r+ic) ;
1644           outval1 += valMask*tex2D(tex_img_inc, j+1, i-r+ic) ;
1645           outval2 += valMask*tex2D(tex_img_inc, j+2, i-r+ic) ;
1646           outval3 += valMask*tex2D(tex_img_inc, j+3, i-r+ic) ;
1647           outval4 += valMask*tex2D(tex_img_inc, j+4, i-r+ic) ;
1648           outval5 += valMask*tex2D(tex_img_inc, j+5, i-r+ic) ;
1649           outval6 += valMask*tex2D(tex_img_inc, j+6, i-r+ic) ;
1650           outval7 += valMask*tex2D(tex_img_inc, j+7, i-r+ic) ;
1651         }
1652
1653   baseIdx = __mul24(i  , j_dim) + j ;
1654   output[ baseIdx++ ] = outval0 ;
1655   output[ baseIdx++ ] = outval1 ;
1656   output[ baseIdx++ ] = outval2 ;
1657   output[ baseIdx++ ] = outval3 ;
1658   output[ baseIdx++ ] = outval4 ;
1659   output[ baseIdx++ ] = outval5 ;
1660   output[ baseIdx++ ] = outval6 ;
1661   output[ baseIdx++ ] = outval7 ;
1662 }
1663
1664
1665 //8 pixels en ligne par thread
1666 __global__ void kernel_convoSepShx8pV(unsigned char *output, int i_dim, int j_dim, int r)
1667 {
1668   int idb, ic, jc, p;
1669   int L = 2*r+1 ;
1670   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
1671   int bdimX = blockDim.x<<3 ;
1672   int tidX = threadIdx.x<<3 ;
1673
1674   // coordonnees absolues du point de base
1675   int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<3 ;
1676   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1677
1678   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
1679
1680
1681   // chargement en smem
1682   int idrow = threadIdx.y*bdimX ;
1683
1684   extern __shared__ unsigned char roi8p[];
1685
1686   // bloc 0 (en haut)
1687   for (p=0; p<8; p++)
1688         roi8p[ idrow  + tidX +p ] = tex2D(tex_img_inc, j+p  , i-r) ;
1689
1690
1691   // bloc 2 ( en bas)
1692   if ( threadIdx.y < L-1 )
1693         {
1694           idrow = (threadIdx.y+blockDim.y)*bdimX ;
1695           for (int p=0; p<8; p++)
1696                 roi8p[ idrow + tidX +p  ] = tex2D( tex_img_inc, j+p  , i+blockDim.y-r ) ;
1697
1698         }
1699   __syncthreads();
1700
1701   // calculs de convolution
1702   // passe verticale
1703   for (ic=0 ; ic<L ; ic++)
1704           {
1705                 int baseRoi = __umul24(ic+threadIdx.y,bdimX) + tidX ;
1706                 float valMask = masque[ ic ] ;
1707                 outval0 += valMask*roi8p[ baseRoi    ] ;
1708                 outval1 += valMask*roi8p[ baseRoi +1 ] ;
1709                 outval2 += valMask*roi8p[ baseRoi +2 ] ;
1710                 outval3 += valMask*roi8p[ baseRoi +3 ] ;
1711                 outval4 += valMask*roi8p[ baseRoi +4 ] ;
1712                 outval5 += valMask*roi8p[ baseRoi +5 ] ;
1713                 outval6 += valMask*roi8p[ baseRoi +6 ] ;
1714                 outval7 += valMask*roi8p[ baseRoi +7 ] ;
1715           }
1716
1717   // 8 pixel par thread --> global mem
1718   output[ idx   ] = outval0 ;
1719   output[ idx+1 ] = outval1 ;
1720   output[ idx+2 ] = outval2 ;
1721   output[ idx+3 ] = outval3 ;
1722   output[ idx+4 ] = outval4 ;
1723   output[ idx+5 ] = outval5 ;
1724   output[ idx+6 ] = outval6 ;
1725   output[ idx+7 ] = outval7 ;
1726 }
1727
1728 __global__ void kernel_convoSepShx8pH(unsigned char *output, int i_dim, int j_dim, int r)
1729 {
1730   int idb, ic, jc, p;
1731   int L = 2*r+1 ;
1732   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
1733   int bdimX = blockDim.x<<3 ;
1734   int tidX = threadIdx.x<<3 ;
1735
1736   // coordonnees absolues du point de base
1737   int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<3 ;
1738   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1739   int j0= __umul24(blockIdx.x,blockDim.x)<<3 ;
1740   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
1741
1742
1743   // chargement en smem
1744   int idrow = threadIdx.y*(bdimX+L-1) ;
1745
1746   extern __shared__ unsigned char roi8p[];
1747
1748   // bloc 0 (a gauche)
1749   for (p=0; p<8; p++)
1750         roi8p[  idrow  + tidX +p ] = tex2D(tex_img_inc, j-r+p  , i) ;
1751
1752   // a droite
1753   if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
1754         {
1755           roi8p[  idrow + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i ) ;
1756           roi8p[  idrow + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i ) ;
1757         }
1758
1759   __syncthreads();
1760
1761   // calculs de convolution
1762   // passe horizontale
1763   for (jc=0 ; jc<L ; jc++)
1764           {
1765                 int baseRoi = idrow + tidX +jc ;
1766                 float valMask = masque[ jc ] ;
1767                 outval0 += valMask*roi8p[  baseRoi    ] ;
1768                 outval1 += valMask*roi8p[  baseRoi +1 ] ;
1769                 outval2 += valMask*roi8p[  baseRoi +2 ] ;
1770                 outval3 += valMask*roi8p[  baseRoi +3 ] ;
1771                 outval4 += valMask*roi8p[  baseRoi +4 ] ;
1772                 outval5 += valMask*roi8p[  baseRoi +5 ] ;
1773                 outval6 += valMask*roi8p[  baseRoi +6 ] ;
1774                 outval7 += valMask*roi8p[  baseRoi +7 ] ;
1775           }
1776
1777   // 1 pixel par thread --> global mem
1778   output[ idx   ] = outval0 ;
1779   output[ idx+1 ] = outval1 ;
1780   output[ idx+2 ] = outval2 ;
1781   output[ idx+3 ] = outval3 ;
1782   output[ idx+4 ] = outval4 ;
1783   output[ idx+5 ] = outval5 ;
1784   output[ idx+6 ] = outval6 ;
1785   output[ idx+7 ] = outval7 ;
1786 }
1787
1788
1789 /*************************************************************************************************
1790  ***********************************************************************************************
1791
1792               FIN DES kERNELS de CONVOLUTION
1793
1794  ***********************************************************************************************
1795  *************************************************************************************************/
1796 // kernel de la libjacket
1797 // Exchange trick: Morgan McGuire, ShaderX 2008
1798 #define s2(a,b)            { float tmp = a; a = min(a,b); b = max(tmp,b); }
1799 #define mn3(a,b,c)         s2(a,b); s2(a,c);
1800 #define mx3(a,b,c)         s2(b,c); s2(a,c);
1801
1802 #define mnmx3(a,b,c)       mx3(a,b,c); s2(a,b);                               // 3 exchanges
1803 #define mnmx4(a,b,c,d)     s2(a,b); s2(c,d); s2(a,c); s2(b,d);                // 4 exchanges
1804 #define mnmx5(a,b,c,d,e)   s2(a,b); s2(c,d); mn3(a,c,e); mx3(b,d,e);          // 6 exchanges
1805 #define mnmx6(a,b,c,d,e,f) s2(a,d); s2(b,e); s2(c,f); mn3(a,b,c); mx3(d,e,f); // 7 exchanges
1806
1807 #define IN(X,Y)  ((0 <= (X) && (X) < nx && 0 <= (Y) && (Y) < ny) ? d_in[(Y)*nx+(X)] : 0)
1808
1809 __global__ static void kernel_medjacket(int nx, int ny, unsigned short*d_out, unsigned short*d_in)
1810 {
1811   int x = __mul24(blockIdx.x , blockDim.x) + threadIdx.x;
1812   int y = __mul24(blockIdx.y , blockDim.y) + threadIdx.y;
1813
1814     // pull top six from image
1815     float v[6] = { IN(x-1, y-1), IN(x  , y-1), IN(x+1, y-1),
1816                    IN(x-1, y  ), IN(x  , y  ), IN(x+1, y  ) };
1817
1818     // with each pass, remove min and max values and add new value
1819     mnmx6(v[0], v[1], v[2], v[3], v[4], v[5]);
1820     v[5] = IN(x-1, y+1); // add new contestant
1821     mnmx5(v[1], v[2], v[3], v[4], v[5]);
1822     v[5] = IN(x  , y+1); // add new contestant
1823     mnmx4(v[2], v[3], v[4], v[5]);
1824     v[5] = IN(x+1, y+1); // add last contenstant
1825     mnmx3(v[3], v[4], v[5]);
1826
1827     // pick the middle one
1828     d_out[y*nx + x] = v[4];
1829 }
1830
1831
1832 /***************************************************************
1833  *   fonction de tri de 2 valeurs entieres (min en premier)
1834  ***************************************************************/
1835 __device__ inline void s(int* a, int* b)
1836 {
1837
1838   int tmp ;
1839   if (*a > *b)
1840         {
1841           tmp = *b ;
1842           *b = *a ;
1843           *a = tmp ;
1844           }
1845 }
1846
1847 __device__ inline void s(unsigned short * a, unsigned short* b)
1848 {
1849
1850   unsigned short tmp ;
1851   if (*a > *b)
1852         {
1853           tmp = *b ;
1854           *b = *a ;
1855           *a = tmp ;
1856           }
1857 }
1858
1859 __device__ inline void s(unsigned char * a, unsigned char * b)
1860 {
1861
1862   unsigned short tmp ;
1863   if (*a > *b)
1864         {
1865           tmp = *b ;
1866           *b = *a ;
1867           *a = tmp ;
1868           }
1869 }
1870
1871 /***************************************************************
1872  *   fonction de min-max d'un tableau de n valeurs entieres
1873  ***************************************************************/
1874 __device__ void minmaxN(unsigned short* v, int n)
1875 {
1876   for (int i=1; i< n; i++)
1877         s( v, v+i) ;
1878   for (int i=n-2; i>0; i--)
1879         s( v+i, v+n-1) ;
1880 }
1881
1882 /***************************************************************
1883  *   fonction de tri naif d'un tableau de n valeurs entieres
1884  ***************************************************************/
1885 __device__ void bubTriN(int * v, int n)
1886 {
1887   for (int i=0; i< n-1; i++)
1888         for (int j=i+1; j<n; j++)
1889           s( v+i, v+j) ;
1890 }
1891
1892
1893 /*****************************************************************************************
1894  *    MACROS pour tri min-max utilisant la fonction de tri s(a,b) ci-dessus
1895  *****************************************************************************************/
1896 #define min3(a, b, c) s(a, b); s(a, c);
1897 #define max3(a, b, c) s(b, c); s(a, c);
1898 #define minmax3(a, b, c) max3(a, b, c); s(a, b); // max 3
1899 #define minmax4(a, b, c, d) s(a, b); s(c, d); s(a, c); s(b, d); //max 4
1900 #define minmax5(a, b, c, d, e) s(a, b); s(c, d); min3(a, c, e); max3(b, d, e); //max 6
1901 #define minmax6(a, b, c, d, e, f) s(a,d); s(b, e); s(c, f); min3(a, b, c); max3(d, e, f); //max 7
1902
1903 /***************************************************************
1904  *   fonction de tri de 9 valeurs entieres
1905  ***************************************************************/
1906 #define bubTriReg9(a, b, c, d, e, f, g, h, i) s(a,b);s(c,d);s(e,f);minmax3(g,h,i);s(a,c);s(e,g);s(a,e);s(b,c);s(d,e);s(f,g);s(b,d);s(f,h);s(b,f);s(c,d);s(e,f);minmax3(g,h,i);s(c,e);s(c,g);s(d,e);s(f,g);s(d,f);s(d,h);s(e,f);minmax3(g,h,i);s(e,g);
1907
1908
1909 __global__ void kernel_bubMedian3( unsigned short*output, int i_dim, int j_dim)
1910 {
1911
1912   // coordonnees absolues du point
1913   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
1914   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
1915
1916   /**************************************************************************
1917    *             tri(s)
1918    **************************************************************************/
1919   int a0, a1, a2, a3, a4, a5, a6, a7, a8 ;
1920
1921   /********************************************************************************
1922    * affectation des valeurs
1923    ********************************************************************************/
1924   a0 = tex2D(tex_img_ins, j-1, i-1) ;
1925   a1 = tex2D(tex_img_ins, j  , i-1) ;
1926   a2 = tex2D(tex_img_ins, j+1, i-1) ;
1927   a3 = tex2D(tex_img_ins, j-1, i) ;
1928   a4 = tex2D(tex_img_ins, j  , i) ;
1929   a5 = tex2D(tex_img_ins, j+1, i) ;
1930   a6 = tex2D(tex_img_ins, j-1, i+1) ;
1931   a7 = tex2D(tex_img_ins, j  , i+1) ;
1932   a8 = tex2D(tex_img_ins, j+1, i+1) ;
1933
1934   //tri selection
1935   bubTriReg9(&a0, &a1, &a2, &a3, &a4, &a5, &a6, &a7, &a8);
1936
1937   //median au milieu !
1938   output[ __mul24(i, j_dim) +j ] = a4 ;
1939 }
1940
1941 __global__ void kernel_median3( unsigned short*output, int i_dim, int j_dim)
1942 {
1943
1944   // coordonnees absolues du point
1945   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
1946   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
1947
1948   /**************************************************************************
1949    *             tri(s)
1950    **************************************************************************/
1951   int a0, a1, a2, a3, a4, a5 ;
1952
1953   /********************************************************************************
1954    * les six premieres valeurs (suffisant pour median 3x3 par forgetfull selection)
1955    ********************************************************************************/
1956   a0 = tex2D(tex_img_ins, j-1, i-1) ;
1957   a1 = tex2D(tex_img_ins, j, i-1) ;
1958   a2 = tex2D(tex_img_ins, j+1, i-1) ;
1959   a3 = tex2D(tex_img_ins, j-1, i) ;
1960   a4 = tex2D(tex_img_ins, j, i) ;
1961   a5 = tex2D(tex_img_ins, j+1, i) ;
1962
1963
1964   //min max aux extremites
1965   minmax6(&a0, &a1, &a2, &a3, &a4, &a5) ;
1966
1967   /********************************************
1968    * les deux valeurs suivantes aux extremites
1969    ********************************************/
1970   a5 = tex2D(tex_img_in, j-1, i+1) ;
1971
1972   minmax5(&a1, &a2, &a3, &a4, &a5) ;
1973
1974   /********************************************
1975    * la derniere valeur a la fin
1976    ********************************************/
1977
1978   a5 = tex2D(tex_img_ins, j, i+1) ;
1979
1980   minmax4(&a2, &a3, &a4, &a5) ;
1981
1982   a5 = tex2D(tex_img_ins, j+1, i+1) ;
1983   minmax3(&a3, &a4, &a5) ;
1984
1985
1986   //median au milieu !
1987   output[ __mul24(i, j_dim) +j ] = a4 ;
1988
1989 }
1990
1991 __global__ void kernel_median3Sh( unsigned short*output, int i_dim, int j_dim)
1992 {
1993   int idb, ic, jc ;
1994   int L = 3 ;
1995
1996   // coordonnees absolues du point de base
1997   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
1998   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1999   int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
2000
2001   // chargement en smem
2002   int idrow = threadIdx.y*(blockDim.x+2) ;
2003
2004   extern __shared__ int medRoi[];
2005
2006   // bloc 0 (en haut à gauche)
2007   medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-1, i-1) ;
2008   // bloc 1 (en haut à droite)...
2009   if ( threadIdx.x < L-1 ) //...ou plutot ce qu'il en manque
2010         medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-1, i-1 ) ;
2011   // bloc 2 ( en bas à gauche)
2012   if ( threadIdx.y < L-1 )
2013         {
2014           idrow = (threadIdx.y+blockDim.y)*(blockDim.x+L-1) ;
2015           medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-1, i+blockDim.y-1 ) ;
2016           //bloc 4 ( en bas à doite )...
2017           if ( threadIdx.x < L-1 ) //...ou ce qu'il en manque
2018                 medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-1, i+blockDim.y-1 ) ;
2019         }
2020   __syncthreads();
2021
2022   /**************************************************************************
2023    *             tri(s)
2024    **************************************************************************/
2025   int a0, a1, a2, a3, a4, a5 ;
2026
2027   /********************************************************************************
2028    * les six premieres valeurs (suffisant pour median 3x3 par forgetfull selection)
2029    ********************************************************************************/
2030   // l'index dans medRoi est du type (threadIdx.y+ic)*(blockDim.X+L-1)+ (threadIdx.x+jc)
2031   // ou ic,jc sont les coordonnees du point dans le noyau 3x3
2032   a0 = medRoi[ (threadIdx.y)*(blockDim.x+L-1)+ (threadIdx.x+1) ] ;
2033   a1 = medRoi[ (threadIdx.y)*(blockDim.x+L-1)+ (threadIdx.x+2)  ] ;
2034   a2 = medRoi[ (threadIdx.y+1)*(blockDim.x+L-1)+ (threadIdx.x+1)] ;
2035   a3 = medRoi[ (threadIdx.y+1)*(blockDim.x+L-1)+ (threadIdx.x+2)] ;
2036   a4 = medRoi[ (threadIdx.y+2)*(blockDim.x+L-1)+ (threadIdx.x+1)] ;
2037   a5 = medRoi[ (threadIdx.y+2)*(blockDim.x+L-1)+ (threadIdx.x+2)] ;
2038
2039
2040   //min max aux extremites
2041   minmax6(&a0, &a1, &a2, &a3, &a4, &a5) ;
2042
2043   /********************************************
2044    * les deux valeurs suivantes aux extremites
2045    ********************************************/
2046   a5 = medRoi[ (threadIdx.y)*(blockDim.x+L-1)+ (threadIdx.x)] ;
2047
2048   minmax5(&a1, &a2, &a3, &a4, &a5) ;
2049
2050   /********************************************
2051    * la derniere valeur a la fin
2052    ********************************************/
2053
2054   a5 = medRoi[ (threadIdx.y+1)*(blockDim.x+L-1)+ (threadIdx.x)] ;
2055
2056   minmax4(&a2, &a3, &a4, &a5) ;
2057
2058   a5 = medRoi[ (threadIdx.y+2)*(blockDim.x+L-1)+ (threadIdx.x)] ;
2059   minmax3(&a3, &a4, &a5) ;
2060
2061
2062   //median au milieu !
2063   output[ __mul24(i, j_dim) +j ] = a4 ;
2064
2065 }
2066
2067 /************************************************
2068  * version qui gere 2 pxels par thread
2069  * Objectif plus déebit en cachant la latence
2070  * minimiser les acces en texture
2071  * defaut = plus de registres
2072  ************************************************/
2073 __global__ void kernel_median3_2pix( unsigned short*output, int i_dim, int j_dim)
2074 {
2075
2076   // coordonnees absolues du point
2077   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
2078   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2079
2080   /**************************************************************************
2081    *             tri(s)
2082    **************************************************************************/
2083   int a0, a1, a2, a3, a4, a5 ;
2084   int b0, b1, b2, b3, b4, b5 ;
2085
2086   /********************************************************************************
2087    * les six premieres valeurs (suffisant pour median 3x3 par forgetfull selection)
2088    ********************************************************************************/
2089   a0 = tex2D(tex_img_ins, j  , i-1) ;
2090   a1 = tex2D(tex_img_ins, j+1, i-1) ;
2091   a2 = tex2D(tex_img_ins, j  , i  ) ;
2092   a3 = tex2D(tex_img_ins, j+1, i  ) ;
2093   a4 = tex2D(tex_img_ins, j  , i+1) ;
2094   a5 = tex2D(tex_img_ins, j+1, i+1) ;
2095
2096
2097   //min max aux extremites
2098   minmax6(&a0, &a1, &a2, &a3, &a4, &a5) ;
2099
2100   b0=a0; b1=a1; b2=a2; b3=a3; b4=a4; b5=a5;
2101
2102   /********************************************
2103    * la valeur suivante au bout
2104    ********************************************/
2105   a5 = tex2D(tex_img_ins, j-1, i  ) ;
2106   b5 = tex2D(tex_img_ins, j+2, i  ) ;
2107
2108   minmax5(&a1, &a2, &a3, &a4, &a5) ;
2109   minmax5(&b1, &b2, &b3, &b4, &b5) ;
2110
2111   /********************************************
2112    * la derniere valeur a la fin
2113    ********************************************/
2114   a5 = tex2D(tex_img_ins, j-1, i-1) ;
2115   b5 = tex2D(tex_img_ins, j+2, i-1) ;
2116
2117   minmax4(&a2, &a3, &a4, &a5) ;
2118   minmax4(&b2, &b3, &b4, &b5) ;
2119
2120   a5 = tex2D(tex_img_ins, j-1, i+1) ;
2121   b5 = tex2D(tex_img_ins, j+2, i+1) ;
2122
2123   minmax3(&a3, &a4, &a5) ;
2124   minmax3(&b3, &b4, &b5) ;
2125
2126   //median au milieu !
2127   output[ __mul24(i, j_dim) +j ] = a4 ;
2128   output[ __mul24(i, j_dim) +j+1 ] = b4 ;
2129
2130 }
2131
2132 __global__ void kernel_median3_2pix8( unsigned char  *output, int i_dim, int j_dim)
2133 {
2134
2135   // coordonnees absolues du point
2136   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
2137   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2138
2139   /**************************************************************************
2140    *             tri(s)
2141    **************************************************************************/
2142   int a0, a1, a2, a3, a4, a5 ;
2143   int b0, b1, b2, b3, b4, b5 ;
2144
2145   /********************************************************************************
2146    * les six premieres valeurs (suffisant pour median 3x3 par forgetfull selection)
2147    ********************************************************************************/
2148   a0 = tex2D(tex_img_inc, j  , i-1) ;
2149   a1 = tex2D(tex_img_inc, j+1, i-1) ;
2150   a2 = tex2D(tex_img_inc, j  , i  ) ;
2151   a3 = tex2D(tex_img_inc, j+1, i  ) ;
2152   a4 = tex2D(tex_img_inc, j  , i+1) ;
2153   a5 = tex2D(tex_img_inc, j+1, i+1) ;
2154
2155
2156   //min max aux extremites
2157   minmax6(&a0, &a1, &a2, &a3, &a4, &a5) ;
2158
2159   b0=a0; b1=a1; b2=a2; b3=a3; b4=a4; b5=a5;
2160
2161   /********************************************
2162    * la valeur suivante au bout
2163    ********************************************/
2164   a5 = tex2D(tex_img_inc, j-1, i  ) ;
2165   b5 = tex2D(tex_img_inc, j+2, i  ) ;
2166
2167   minmax5(&a1, &a2, &a3, &a4, &a5) ;
2168   minmax5(&b1, &b2, &b3, &b4, &b5) ;
2169
2170   /********************************************
2171    * la derniere valeur a la fin
2172    ********************************************/
2173   a5 = tex2D(tex_img_inc, j-1, i-1) ;
2174   b5 = tex2D(tex_img_inc, j+2, i-1) ;
2175
2176   minmax4(&a2, &a3, &a4, &a5) ;
2177   minmax4(&b2, &b3, &b4, &b5) ;
2178
2179   a5 = tex2D(tex_img_inc, j-1, i+1) ;
2180   b5 = tex2D(tex_img_inc, j+2, i+1) ;
2181
2182   minmax3(&a3, &a4, &a5) ;
2183   minmax3(&b3, &b4, &b5) ;
2184
2185   //median au milieu !
2186   output[ __mul24(i, j_dim) +j ] = a4 ;
2187   output[ __mul24(i, j_dim) +j+1 ] = b4 ;
2188
2189 }
2190
2191
2192 #define minmax7(a,b,c,d,e,f,g) s(a,b);s(c,d);s(e,f);s(a,c);s(e,g);s(a,e);s(b,d);s(f,g);s(d,g);//max 9
2193 #define minmax8(a,b,c,d,e,f,g,h) s(a,b);s(c,d);s(e,f);s(g,h);s(a,c);s(e,g);s(a,e);s(b,d);s(f,h);s(d,h);//max 10
2194 #define minmax9(a,b,c,d,e,f,g,h,i) s(a,b);s(c,d);s(e,f);minmax3(g,h,i);s(a,c);s(e,g);s(a,e);s(b,d);s(f,i);s(d,i);//max 12
2195 #define minmax10(a,b,c,d,e,f,g,h,i,j) s(a,b);s(c,d);s(e,f);s(g,h);s(i,j);s(a,c);min3(e,g,i);s(a,e);s(b,d);max3(f,h,j);s(d,j);//max 13
2196 #define minmax11(a,b,c,d,e,f,g,h,i,j,k) s(a,b);s(c,d);s(e,f);s(g,h);min3(i,j,k);s(a,c);min3(e,g,i);s(a,e);s(b,d);max3(f,h,k);s(d,k);//max 14
2197 #define minmax12(a,b,c,d,e,f,g,h,i,j,k,l) s(a,b);s(c,d);s(e,f);s(g,h);s(i,j);s(k,l);s(a,c);s(e,g);s(i,k);min3(a,e,i);s(b,d);s(f,h);s(j,l);max3(d,h,l); //max 16
2198 #define minmax13(a,b,c,d,e,f,g,h,i,j,k,l,m) s(a,b);s(c,d);s(e,f);s(g,h);s(i,j);min3(k,l,m);s(a,c);s(e,g);s(i,k);min3(a,e,i);s(b,d);s(f,h);max3(j,l,m);max3(d,h,m); // max 18 !!
2199 #define min14(a,b,c,d,e,f,g,h,i,j,k,l,m,n) s(a,b);s(c,d);s(e,f);s(g,h);s(i,j);s(k,l);s(m,n);s(a,c);s(e,g);s(i,k);s(a,e);s(i,m);s(a,i); // 13 mouvts max
2200 #define max14(a,b,c,d,e,f,g,h,i,j,k,l,m,n) s(a,b);s(c,d);s(e,f);s(g,h);s(i,j);s(k,l);s(m,n);s(b,d);s(f,h);s(j,l);s(d,h);s(l,n);s(h,n); // 13 mouvts max
2201 #define minmax14(a,b,c,d,e,f,g,h,i,j,k,l,m,n) min14(a,b,c,d,e,f,g,h,i,j,k,l,m,n); s(b,d);s(f,h);s(j,l);s(d,h);s(l,n);s(h,n);           // 19 mouvts max
2202
2203 #define m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o) s(a,b);s(c,d);s(e,f);s(g,h);s(i,j);s(k,l);s(a,c);s(e,g);s(i,k);s(b,d);s(f,h);s(j,l); // max 12
2204 #define minmax15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);minmax3(m,n,o);s(a,e);s(i,m);s(d,h);s(l,o);s(a,i);s(h,o); //max 21
2205 #define minmax16(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(m,o);s(a,e);s(i,m);s(n,p);s(d,h);s(l,p);s(a,i);s(h,p);//max 22
2206 #define minmax17(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);minmax3(o,p,q);s(a,e);s(d,h);s(i,m);min3(a,i,o);max3(h,n,q);//max 24
2207 #define minmax18(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);minmax3(m,o,q);minmax3(n,p,r);s(a,e);s(m,n);s(q,r);s(i,m);s(d,h);s(l,r);s(a,i);s(h,r); //max 29 !!!
2208 #define minmax19(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);minmax3(q,r,ss);s(m,o);s(n,p);s(a,e);s(i,m);s(d,h);s(l,p);min3(a,i,q);max3(h,p,ss);//max 28
2209 #define minmax20(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss,t) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);s(ss,t);s(m,o);s(q,ss);s(n,p);s(r,t);s(a,e);s(i,m);s(d,h);s(l,p);min3(a,i,q);max3(h,p,t);//max 29
2210 #define minmax21(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss,t,u) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);minmax3(ss,t,u);s(m,o);s(q,ss);s(n,p);s(r,u);s(a,e);s(i,m);s(d,h);s(l,p);min3(a,i,q);max3(h,p,u); //max 30
2211 #define minmax22(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss,t,u,v) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);s(ss,t);s(u,v);s(m,o);min3(q,ss,u);s(n,p);max3(r,t,v);s(a,e);s(i,m);s(d,h);s(l,p);min3(a,i,q);max3(h,p,v); //max31
2212 #define minmax23(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss,t,u,v,w) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);s(ss,t);minmax3(u,v,w);s(m,o);s(q,ss);s(n,p);s(r,t);s(a,e);s(i,m);s(q,u);s(d,h);s(l,p);s(t,w);min3(a,i,q);max3(h,p,w); //max 33
2213 #define minmax24(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss,t,u,v,w,x) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);s(ss,t);s(u,v);s(w,x);s(m,o);s(q,ss);s(u,w);s(n,p);s(r,t);s(v,x);s(a,e);s(i,m);s(q,u);s(d,h);s(l,p);s(t,x);min3(a,i,q);max3(h,p,x); //max 34
2214 #define minmax25(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss,t,u,v,w,x,y) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);s(ss,t);s(u,v);minmax3(w,x,y);s(m,o);s(q,ss);s(u,w);s(n,p);s(r,t);s(v,y);s(a,e);s(i,m);s(q,u);s(d,h);s(l,p);s(t,y);min3(a,i,q);max3(h,p,y); //max 36
2215 #define minmax26(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss,t,u,v,w,x,y,z) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);s(ss,t);s(u,v);s(w,x);s(y,z);s(m,o);s(q,ss);min3(u,w,y);s(n,p);s(r,t);max3(v,x,z);s(a,e);s(i,m);s(q,u);s(d,h);s(l,p);s(t,z);min3(a,i,q);max3(h,p,z);//max 37
2216
2217
2218 __global__ void kernel_median5( unsigned short*output, int i_dim, int j_dim)
2219 {
2220
2221   // coordonnees absolues du point
2222   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
2223   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2224
2225   /**************************************************************************
2226    *             tri(s)
2227    **************************************************************************/
2228   int a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13 ;
2229
2230   /********************************************************************************
2231    * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
2232    ********************************************************************************/
2233   //premiere ligne
2234   a0  = tex2D(tex_img_ins, j-2, i-2) ;
2235   a1  = tex2D(tex_img_ins, j-1, i-2) ;
2236   a2  = tex2D(tex_img_ins, j  , i-2) ;
2237   a3  = tex2D(tex_img_ins, j+1, i-2) ;
2238   a4  = tex2D(tex_img_ins, j+2, i-2) ;
2239   //deuxieme ligne
2240   a5  = tex2D(tex_img_ins, j-2, i-1) ;
2241   a6  = tex2D(tex_img_ins, j-1, i-1) ;
2242   a7  = tex2D(tex_img_ins, j  , i-1) ;
2243   a8  = tex2D(tex_img_ins, j+1, i-1) ;
2244   a9  = tex2D(tex_img_ins, j+2, i-1) ;
2245   //les 4 premiers de la troisieme ligne
2246   a10 = tex2D(tex_img_ins, j-2, i) ;
2247   a11 = tex2D(tex_img_ins, j-1, i) ;
2248   a12 = tex2D(tex_img_ins, j  , i) ;
2249   a13 = tex2D(tex_img_ins, j+1, i) ;
2250
2251   //min max aux extremites
2252   minmax14(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2253
2254   //chargement valeurs suivante (15)
2255   a13 = tex2D(tex_img_ins, j+2, i);
2256   //minmax aux extremites
2257   minmax13(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2258
2259
2260   //chargement valeur suivante (16)
2261   a13 = tex2D(tex_img_ins, j-2, i+1);
2262   //minmax aux extremites
2263   minmax12(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2264
2265   //chargement valeur suivante (17)
2266   a13 = tex2D(tex_img_ins, j-1, i+1);
2267   //minmax aux extremites
2268   minmax11(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2269
2270   //chargement valeur suivante (18)
2271   a13 = tex2D(tex_img_ins, j  , i+1);
2272   //minmax aux extremites
2273   minmax10(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2274
2275   //chargement valeur suivante (19)
2276   a13 = tex2D(tex_img_ins, j+1, i+1);
2277   //minmax aux extremites
2278   minmax9(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2279
2280   //chargement valeur suivante (20)
2281   a13 = tex2D(tex_img_ins, j+2, i+1);
2282   //minmax aux extremites
2283   minmax8(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2284
2285   //chargement valeur suivante (21)
2286   a13 = tex2D(tex_img_ins, j-2, i+2);
2287   //minmax aux extremites
2288   minmax7(&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2289
2290   //chargement valeur suivante (22)
2291   a13 = tex2D(tex_img_ins, j-1, i+2);
2292   //minmax aux extremites
2293   minmax6(&a8,&a9,&a10,&a11,&a12,&a13);
2294
2295   //chargement valeur suivante (23)
2296   a13 = tex2D(tex_img_ins, j  , i+2);
2297   //minmax aux extremites
2298   minmax5(&a9,&a10,&a11,&a12,&a13);
2299
2300   //chargement valeur suivante (24)
2301   a13 = tex2D(tex_img_ins, j+1  , i+2);
2302   //minmax aux extremites
2303   minmax4(&a10,&a11,&a12,&a13);
2304
2305   //chargement valeur suivante (25)
2306   a13 = tex2D(tex_img_ins, j+2  , i+2);
2307   //minmax aux extremites
2308   minmax3(&a11,&a12,&a13);
2309
2310
2311   //median au milieu !
2312   output[ __mul24(i, j_dim) +j ] = a12 ;
2313
2314 }
2315
2316 __global__ void kernel_median5_sh( unsigned short*output, int i_dim, int j_dim)
2317 {
2318
2319   // coordonnees absolues du point
2320   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
2321   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2322   int L = 5 ;
2323   int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
2324
2325   // chargement en smem
2326   int idrow = threadIdx.y*(blockDim.x+4) ;
2327
2328   extern __shared__ int medRoi[];
2329
2330   // bloc 0 (en haut à gauche)
2331   medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-2, i-2) ;
2332   // bloc 1 (en haut à droite)...
2333   if ( threadIdx.x < L-1 ) //...ou plutot ce qu'il en manque
2334         medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-2, i-2 ) ;
2335   // bloc 2 ( en bas à gauche)
2336   if ( threadIdx.y < L-1 )
2337         {
2338           idrow = (threadIdx.y+blockDim.y)*(blockDim.x+L-1) ;
2339           medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-2, i+blockDim.y-2 ) ;
2340           //bloc 4 ( en bas à doite )...
2341           if ( threadIdx.x < L-1 ) //...ou ce qu'il en manque
2342                 medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-2, i+blockDim.y-2 ) ;
2343         }
2344   __syncthreads();
2345
2346   /**************************************************************************
2347    *             tri(s)
2348    **************************************************************************/
2349   int a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13 ;
2350
2351   /********************************************************************************
2352    * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
2353    ********************************************************************************/
2354   // l'index dans medRoi est du type (threadIdx.y+ic)*(blockDim.x+L-1)+ (threadIdx.x+jc)
2355   // ou ic,jc sont les coordonnees du point dans le noya
2356   //premiere ligne
2357   a0  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i-2) ;
2358   a1  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i-2) ;
2359   a2  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i-2) ;
2360   a3  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i-2) ;
2361   a4  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i-2) ;
2362   //deuxieme ligne
2363   a5  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i-1) ;
2364   a6  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i-1) ;
2365   a7  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i-1) ;
2366   a8  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i-1) ;
2367   a9  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i-1) ;
2368   //les 4 premiers de la troisieme ligne
2369   a10 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i) ;
2370   a11 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i) ;
2371   a12 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i) ;
2372   a13 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i) ;
2373
2374   //min max aux extremites
2375   minmax14(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2376
2377   //chargement valeurs suivante (15)
2378   a13 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i);
2379   //minmax aux extremites
2380   minmax13(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2381
2382
2383   //chargement valeur suivante (16)
2384   a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i+1);
2385   //minmax aux extremites
2386   minmax12(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2387
2388   //chargement valeur suivante (17)
2389   a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i+1);
2390   //minmax aux extremites
2391   minmax11(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2392
2393   //chargement valeur suivante (18)
2394   a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i+1);
2395   //minmax aux extremites
2396   minmax10(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2397
2398   //chargement valeur suivante (19)
2399   a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i+1);
2400   //minmax aux extremites
2401   minmax9(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2402
2403   //chargement valeur suivante (20)
2404   a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i+1);
2405   //minmax aux extremites
2406   minmax8(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2407
2408   //chargement valeur suivante (21)
2409   a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i+2);
2410   //minmax aux extremites
2411   minmax7(&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2412
2413   //chargement valeur suivante (22)
2414   a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i+2);
2415   //minmax aux extremites
2416   minmax6(&a8,&a9,&a10,&a11,&a12,&a13);
2417
2418   //chargement valeur suivante (23)
2419   a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i+2);
2420   //minmax aux extremites
2421   minmax5(&a9,&a10,&a11,&a12,&a13);
2422
2423   //chargement valeur suivante (24)
2424   a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1  , i+2);
2425   //minmax aux extremites
2426   minmax4(&a10,&a11,&a12,&a13);
2427
2428   //chargement valeur suivante (25)
2429   a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2  , i+2);
2430   //minmax aux extremites
2431   minmax3(&a11,&a12,&a13);
2432
2433
2434   //median au milieu !
2435   output[ __mul24(i, j_dim) +j ] = a12 ;
2436
2437 }
2438
2439
2440 __global__ void kernel_median5_sh_loc( unsigned short*output, int i_dim, int j_dim)
2441 {
2442
2443   // coordonnees absolues du point
2444   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
2445   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2446   // int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
2447
2448   // chargement en smem
2449   int idrow = threadIdx.y*(blockDim.x+4) ;
2450
2451   extern __shared__ int medRoi[];
2452
2453   // bloc 0 (en haut à gauche)
2454   medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-2, i-2) ;
2455   // bloc 1 (en haut à droite)...
2456   if ( threadIdx.x < 4 ) //...ou plutot ce qu'il en manque
2457         medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-2, i-2 ) ;
2458   // bloc 2 ( en bas à gauche)
2459   if ( threadIdx.y < 4 )
2460         {
2461           idrow = (threadIdx.y+blockDim.y)*(blockDim.x+4) ;
2462           medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-2, i+blockDim.y-2 ) ;
2463           //bloc 4 ( en bas à doite )...
2464           if ( threadIdx.x < 4 ) //...ou ce qu'il en manque
2465                 medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-2, i+blockDim.y-2 ) ;
2466         }
2467   __syncthreads();
2468
2469   /**************************************************************************
2470    *             tri(s)
2471    **************************************************************************/
2472   /********************************************************************************
2473    * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
2474    ********************************************************************************/
2475   // l'index dans medRoi est du type (threadIdx.y+ic)*(blockDim.x+L-1)+ (threadIdx.x+jc)
2476   // ou ic,jc sont les coordonnees du point dans le noyau
2477
2478   int a[] = { medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x ) ],
2479                           medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x +1) ],
2480                           medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x +2) ],
2481                           medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x +3) ],
2482                           medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x +4) ],
2483                           medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x ) ],
2484                           medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x +1) ],
2485                           medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x +2) ],
2486                           medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x +3) ],
2487                           medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x +4) ],
2488                           medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x ) ],
2489                           medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x +1) ],
2490                           medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x +2) ],
2491                           medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x +3) ]};
2492
2493
2494   //min max aux extremites
2495   minmax14(&a[0],&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2496
2497   //chargement valeurs suivante (15)
2498   a[13] = medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i);
2499   //minmax aux extremites
2500   minmax13(&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2501
2502
2503   //chargement valeur suivante (16)
2504   a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i+1);
2505   //minmax aux extremites
2506   minmax12(&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2507
2508   //chargement valeur suivante (17)
2509   a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i+1);
2510   //minmax aux extremites
2511   minmax11(&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2512
2513   //chargement valeur suivante (18)
2514   a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i+1);
2515   //minmax aux extremites
2516   minmax10(&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2517
2518   //chargement valeur suivante (19)
2519   a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i+1);
2520   //minmax aux extremites
2521   minmax9(&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2522
2523   //chargement valeur suivante (20)
2524   a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i+1);
2525   //minmax aux extremites
2526   minmax8(&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2527
2528   //chargement valeur suivante (21)
2529   a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i+2);
2530   //minmax aux extremites
2531   minmax7(&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2532
2533   //chargement valeur suivante (22)
2534   a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i+2);
2535   //minmax aux extremites
2536   minmax6(&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2537
2538   //chargement valeur suivante (23)
2539   a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i+2);
2540   //minmax aux extremites
2541   minmax5(&a[9],&a[10],&a[11],&a[12],&a[13]);
2542
2543   //chargement valeur suivante (24)
2544   a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1  , i+2);
2545   //minmax aux extremites
2546   minmax4(&a[10],&a[11],&a[12],&a[13]);
2547
2548   //chargement valeur suivante (25)
2549   a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2  , i+2);
2550   //minmax aux extremites
2551   minmax3(&a[11],&a[12],&a[13]);
2552
2553
2554   //median au milieu !
2555   output[ __mul24(i, j_dim) +j ] = a[12] ;
2556
2557 }
2558
2559 __global__ void kernel_median5_sh_loc_2p( unsigned short*output, int i_dim, int j_dim)
2560 {
2561
2562   // coordonnees absolues du point
2563   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
2564   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2565   // int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
2566   int bdimX = blockDim.x<<1 ;
2567   int tidX = threadIdx.x<<1 ;
2568
2569   // chargement en smem
2570   int idrow = threadIdx.y*(bdimX+2*2) ;
2571
2572   extern __shared__ int medRoi[];
2573
2574   // bloc 0 (en haut à gauche)
2575   medRoi[ idrow  + tidX ] = tex2D(tex_img_ins, j-2, i-2) ;
2576   medRoi[ idrow  + tidX +1] = tex2D(tex_img_ins, j-1, i-2) ;
2577   // bloc 1 (en haut à droite)...
2578   if ( tidX < 4 ) //...ou plutot ce qu'il en manque
2579         {
2580           medRoi[ idrow + tidX + bdimX ] = tex2D( tex_img_ins, j+bdimX-2, i-2 ) ;
2581           medRoi[ idrow + tidX + bdimX +1] = tex2D( tex_img_ins, j+bdimX-1, i-2 ) ;
2582         }
2583   // bloc 2 ( en bas à gauche)
2584   if ( threadIdx.y < 4 )
2585         {
2586           idrow = (threadIdx.y+blockDim.y)*(bdimX+4) ;
2587           medRoi[ idrow + tidX ] = tex2D( tex_img_ins, j-2, i+blockDim.y-2 ) ;
2588           medRoi[ idrow + tidX +1] = tex2D( tex_img_ins, j-1, i+blockDim.y-2 ) ;
2589           //bloc 4 ( en bas à doite )...
2590           if ( tidX < 4 ) //...ou ce qu'il en manque
2591                 {
2592                   medRoi[ idrow + tidX + bdimX ] = tex2D( tex_img_ins, j+bdimX-2, i+blockDim.y-2 ) ;
2593                   medRoi[ idrow + tidX + bdimX +1] = tex2D( tex_img_ins, j+bdimX-1, i+blockDim.y-2 ) ;
2594                 }
2595         }
2596   __syncthreads();
2597
2598   /**************************************************************************
2599    *             tri(s)
2600    **************************************************************************/
2601   /********************************************************************************
2602    * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
2603    ********************************************************************************/
2604   // l'index dans medRoi est du type (threadIdx.y+ic)*(bdimX+L-1)+ (tidX+jc)
2605   // ou ic,jc sont les coordonnees du point dans le noyau
2606
2607   //on commence le tri avec les pixels du voisinage commun aux 2 pixels à calculer
2608
2609   int a[] = { medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +1) ],
2610                           medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +2) ],
2611                           medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +3) ],
2612                           medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +4) ],
2613
2614                           medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +1) ],
2615                           medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +2) ],
2616                           medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +3) ],
2617                           medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +4) ],
2618
2619                           medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +1) ],
2620                           medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +2) ],
2621                           medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +3) ],
2622                           medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +4) ],
2623
2624                           medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +1) ],
2625                           medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +2) ]
2626   };
2627
2628   //min max aux extremites
2629   minmax14(&a[0],&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2630
2631   //chargement valeurs suivante (15)
2632   a[13] = medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +3) ] ;
2633   //minmax aux extremites
2634   minmax13(&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2635
2636
2637   //chargement valeur suivante (16)
2638   a[13] = medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +4) ] ;
2639   //minmax aux extremites
2640   minmax12(&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2641
2642   //chargement valeur suivante (17)
2643   a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +1) ] ;
2644   //minmax aux extremites
2645   minmax11(&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2646
2647   //chargement valeur suivante (18)
2648   a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +2) ] ;
2649   //minmax aux extremites
2650   minmax10(&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2651
2652   //chargement valeur suivante (19)
2653   a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +3) ] ;
2654   //minmax aux extremites
2655   minmax9(&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2656
2657   //chargement valeur suivante (20)
2658   a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +4) ] ;
2659   //minmax aux extremites
2660   minmax8(&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2661
2662   // fin des pixels communs. la suite doit se faire en 2 tris distincts.
2663   int b[] = { a[7], a[8], a[9], a[10], a[11], a[12], a[13]};
2664
2665   //chargement valeur suivante (21)
2666   a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX ) ] ;
2667   b[6 ] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +5) ] ;
2668   //minmax aux extremites
2669   minmax7(&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2670   minmax7(&b[0],&b[1],&b[2],&b[3],&b[4],&b[5],&b[6]);
2671
2672   //chargement valeur suivante (22)
2673   a[13] = medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX ) ] ;
2674   b[6 ] = medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +5) ] ;
2675   //minmax aux extremites
2676   minmax6(&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
2677   minmax6(&b[1],&b[2],&b[3],&b[4],&b[5],&b[6]);
2678
2679   //chargement valeur suivante (23)
2680   a[13] = medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX ) ] ;
2681   b[6 ] = medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +5) ] ;
2682   //minmax aux extremites
2683   minmax5(&a[9],&a[10],&a[11],&a[12],&a[13]);
2684   minmax5(&b[2],&b[3],&b[4],&b[5],&b[6]);
2685
2686   //chargement valeur suivante (24)
2687   a[13] = medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX ) ] ;
2688   b[6 ] = medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +5) ] ;
2689
2690   //minmax aux extremites
2691   minmax4(&a[10],&a[11],&a[12],&a[13]);
2692   minmax4(&b[3],&b[4],&b[5],&b[6]);
2693
2694   //chargement valeur suivante (25)
2695   a[13] = medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX ) ] ;
2696   b[6 ] = medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +5) ] ;
2697   //minmax aux extremites
2698   minmax3(&a[11],&a[12],&a[13]);
2699   minmax3(&b[4],&b[5],&b[6]);
2700
2701   //median au milieu !
2702   output[ __mul24(i, j_dim) +j ] = a[12] ;
2703   output[ __mul24(i, j_dim) +j+1 ] = b[5] ;
2704
2705 }
2706
2707
2708 __global__ void kernel_medianR( unsigned short * output, int i_dim, int j_dim, int r)
2709 {
2710   // coordonnees absolues du point
2711   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
2712   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2713   unsigned short cpt ;
2714   unsigned short ic, jc ;
2715
2716   unsigned short hist[256] ;
2717   for (ic =0; ic<256; ic++) hist[ic]=0; // init histogramme
2718   //for (ic =0; ic<25; ic++) hist[ic]=1; // init histogramme
2719
2720   //generation histogramme
2721   for(ic=i-r; ic<=i+r; ic++ )
2722         for(jc=j-r; jc<=j+r; jc++)
2723           hist[tex2D(tex_img_ins, jc, ic)]++ ;
2724
2725   //parcours histogramme
2726   cpt = 0 ;
2727   for(ic=0; ic<256; ic++)
2728         {
2729           cpt += hist[ic] ;
2730           //selection de la valeur du percentile (ici 50%=SUM/2)
2731           if ( cpt > ((2*r+1)*(2*r+1))>>1 ) break ;
2732         }
2733
2734   output[ __mul24(i, j_dim) +j ] =  ic ;
2735 }
2736
2737 __global__ void kernel_medianRSH( unsigned short * output, int i_dim, int j_dim, int r)
2738 {
2739   // coordonnees absolues du point
2740   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
2741   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2742   unsigned short cpt ;
2743   unsigned short ic, jc ;
2744
2745   // chargement en smem
2746   int idrow = threadIdx.y*(blockDim.x+2*r) ;
2747
2748   extern __shared__ int medRoi[];
2749
2750   // bloc 0 (en haut à gauche)
2751   medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-r, i-r) ;
2752   // bloc 1 (en haut à droite)...
2753   if ( threadIdx.x < 2*r ) //...ou plutot ce qu'il en manque
2754         medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-r, i-r ) ;
2755   // bloc 2 ( en bas à gauche)
2756   if ( threadIdx.y < 2*r )
2757         {
2758           idrow = (threadIdx.y+blockDim.y)*(blockDim.x+2*r) ;
2759           medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-r, i+blockDim.y-r ) ;
2760           //bloc 4 ( en bas à doite )...
2761           if ( threadIdx.x < 2*r ) //...ou ce qu'il en manque
2762                 medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-r, i+blockDim.y-r ) ;
2763         }
2764   __syncthreads();
2765
2766   int hist[256] ;
2767   for (ic =0; ic<256; ic++) hist[ic]=0; // init histogramme
2768
2769   //generation histogramme
2770   for(ic=0; ic<=2*r; ic++ )
2771         for(jc=0; jc<=2*r; jc++)
2772           hist[medRoi[(threadIdx.y+ic)*(blockDim.x+2*r)+ (threadIdx.x+jc)]]++ ;
2773
2774   //parcours histogramme
2775   cpt = 0 ;
2776   for(ic=0; ic<256; ic++)
2777         {
2778           cpt += hist[ic] ;
2779           //selection de la valeur du percentile (ici 50%=SUM/2)
2780           if ( cpt > (((2*r+1)*(2*r+1))>>1) ) break ;
2781         }
2782
2783   output[ __mul24(i, j_dim) +j ] =  ic ;
2784 }
2785
2786
2787 __global__ void kernel_median5_2pix( unsigned short*output, int i_dim, int j_dim)
2788 {
2789
2790   // coordonnees absolues du point
2791   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
2792   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2793
2794   /**************************************************************************
2795    *             tri(s)
2796    **************************************************************************/
2797   int a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13 ;
2798   int b7, b8, b9, b10, b11, b12, b13 ;
2799
2800   /********************************************************************************
2801    * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
2802    ********************************************************************************/
2803   //premiere ligne
2804   a0  = tex2D(tex_img_ins, j-1, i-2) ;
2805   a1  = tex2D(tex_img_ins, j  , i-2) ;
2806   a2  = tex2D(tex_img_ins, j+1, i-2) ;
2807   a3  = tex2D(tex_img_ins, j+2, i-2) ;
2808   //deuxieme ligne
2809   a4  = tex2D(tex_img_ins, j-1, i-1) ;
2810   a5  = tex2D(tex_img_ins, j  , i-1) ;
2811   a6  = tex2D(tex_img_ins, j+1, i-1) ;
2812   a7  = tex2D(tex_img_ins, j+2, i-1) ;
2813   //troisieme ligne
2814   a8  = tex2D(tex_img_ins, j-1, i) ;
2815   a9  = tex2D(tex_img_ins, j  , i) ;
2816   a10 = tex2D(tex_img_ins, j+1, i) ;
2817   a11 = tex2D(tex_img_ins, j+2, i) ;
2818   //les 2 premiers de la 4eme ligne
2819   a12 = tex2D(tex_img_ins, j-1, i+1) ;
2820   a13 = tex2D(tex_img_ins, j  , i+1) ;
2821
2822   //min max aux extremites
2823   minmax14(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2824
2825   //chargement valeurs suivante (15)
2826   a13 = tex2D(tex_img_ins, j+1, i+1);
2827   //minmax aux extremites
2828   minmax13(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2829
2830   //chargement valeur suivante (16)
2831   a13 = tex2D(tex_img_ins, j+2, i+1);
2832   //minmax aux extremites
2833   minmax12(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2834
2835   //chargement valeur suivante (17)
2836   a13 = tex2D(tex_img_ins, j-1, i+2);
2837   //minmax aux extremites
2838   minmax11(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2839
2840   //chargement valeur suivante (18)
2841   a13 = tex2D(tex_img_ins, j  , i+2);
2842   //minmax aux extremites
2843   minmax10(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2844
2845   //chargement valeur suivante (19)
2846   a13 = tex2D(tex_img_ins, j+1, i+2);
2847   //minmax aux extremites
2848   minmax9(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2849
2850   //chargement valeur suivante (20)
2851   a13 = tex2D(tex_img_ins, j+2, i+2);
2852   //minmax aux extmites
2853   minmax8(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2854
2855   // fin des pixels voisins communs deux pixels centraux
2856   b7=a7; b8=a8; b9=a9; b10=a10; b11=a11; b12=a12; b13=a13;
2857
2858   //chargement valeur suivante (21)
2859   a13 = tex2D(tex_img_ins, j-2, i-2);
2860   b13 = tex2D(tex_img_ins, j+3, i-2);
2861   //minmax aux extremites
2862   minmax7(&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2863   minmax7(&b7,&b8,&b9,&b10,&b11,&b12,&b13);
2864
2865   //chargement valeur suivante (22)
2866   a13 = tex2D(tex_img_ins, j-2, i-1);
2867   b13 = tex2D(tex_img_ins, j+3, i-1);
2868   //minmax aux extremites
2869   minmax6(&a8,&a9,&a10,&a11,&a12,&a13);
2870   minmax6(&b8,&b9,&b10,&b11,&b12,&b13);
2871
2872   //chargement valeur suivante (23)
2873   a13 = tex2D(tex_img_ins, j-2, i  );
2874   b13 = tex2D(tex_img_ins, j+3, i  );
2875   //minmax aux extremites
2876   minmax5(&a9,&a10,&a11,&a12,&a13);
2877   minmax5(&b9,&b10,&b11,&b12,&b13);
2878
2879   //chargement valeur suivante (24)
2880   a13 = tex2D(tex_img_ins, j-2, i+1);
2881   b13 = tex2D(tex_img_ins, j+3, i+1);
2882   //minmax aux extremites
2883   minmax4(&a10,&a11,&a12,&a13);
2884   minmax4(&b10,&b11,&b12,&b13);
2885
2886   //chargement valeur suivante (25)
2887   a13 = tex2D(tex_img_ins, j-2, i+2);
2888   b13 = tex2D(tex_img_ins, j+3, i+2);
2889   //minmax aux extremites
2890   minmax3(&a11,&a12,&a13);
2891   minmax3(&b11,&b12,&b13);
2892
2893   //median au milieu !
2894   output[ __mul24(i, j_dim) +j ] = a12 ;
2895   output[ __mul24(i, j_dim) +j+1 ] = b12 ;
2896
2897 }
2898
2899
2900 __global__ void kernel_median5_2pix( unsigned char  *output, int i_dim, int j_dim)
2901 {
2902
2903   // coordonnees absolues du point
2904   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
2905   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2906
2907   /**************************************************************************
2908    *             tri(s)
2909    **************************************************************************/
2910   int a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13 ;
2911   int b7, b8, b9, b10, b11, b12, b13 ;
2912
2913   /********************************************************************************
2914    * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
2915    ********************************************************************************/
2916   //premiere ligne
2917   a0  = tex2D(tex_img_inc, j-1, i-2) ;
2918   a1  = tex2D(tex_img_inc, j  , i-2) ;
2919   a2  = tex2D(tex_img_inc, j+1, i-2) ;
2920   a3  = tex2D(tex_img_inc, j+2, i-2) ;
2921   //deuxieme ligne
2922   a4  = tex2D(tex_img_inc, j-1, i-1) ;
2923   a5  = tex2D(tex_img_inc, j  , i-1) ;
2924   a6  = tex2D(tex_img_inc, j+1, i-1) ;
2925   a7  = tex2D(tex_img_inc, j+2, i-1) ;
2926   //troisieme ligne
2927   a8  = tex2D(tex_img_inc, j-1, i) ;
2928   a9  = tex2D(tex_img_inc, j  , i) ;
2929   a10 = tex2D(tex_img_inc, j+1, i) ;
2930   a11 = tex2D(tex_img_inc, j+2, i) ;
2931   //les 2 premiers de la 4eme ligne
2932   a12 = tex2D(tex_img_inc, j-1, i+1) ;
2933   a13 = tex2D(tex_img_inc, j  , i+1) ;
2934
2935   //min max aux extremites
2936   minmax14(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2937
2938   //chargement valeurs suivante (15)
2939   a13 = tex2D(tex_img_inc, j+1, i+1);
2940   //minmax aux extremites
2941   minmax13(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2942
2943   //chargement valeur suivante (16)
2944   a13 = tex2D(tex_img_inc, j+2, i+1);
2945   //minmax aux extremites
2946   minmax12(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2947
2948   //chargement valeur suivante (17)
2949   a13 = tex2D(tex_img_inc, j-1, i+2);
2950   //minmax aux extremites
2951   minmax11(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2952
2953   //chargement valeur suivante (18)
2954   a13 = tex2D(tex_img_inc, j  , i+2);
2955   //minmax aux extremites
2956   minmax10(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2957
2958   //chargement valeur suivante (19)
2959   a13 = tex2D(tex_img_inc, j+1, i+2);
2960   //minmax aux extremites
2961   minmax9(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2962
2963   //chargement valeur suivante (20)
2964   a13 = tex2D(tex_img_inc, j+2, i+2);
2965   //minmax aux extmites
2966   minmax8(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2967
2968   // fin des pixels voisins communs deux pixels centraux
2969   b7=a7; b8=a8; b9=a9; b10=a10; b11=a11; b12=a12; b13=a13;
2970
2971   //chargement valeur suivante (21)
2972   a13 = tex2D(tex_img_inc, j-2, i-2);
2973   b13 = tex2D(tex_img_inc, j+3, i-2);
2974   //minmax aux extremites
2975   minmax7(&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2976   minmax7(&b7,&b8,&b9,&b10,&b11,&b12,&b13);
2977
2978   //chargement valeur suivante (22)
2979   a13 = tex2D(tex_img_inc, j-2, i-1);
2980   b13 = tex2D(tex_img_inc, j+3, i-1);
2981   //minmax aux extremites
2982   minmax6(&a8,&a9,&a10,&a11,&a12,&a13);
2983   minmax6(&b8,&b9,&b10,&b11,&b12,&b13);
2984
2985   //chargement valeur suivante (23)
2986   a13 = tex2D(tex_img_inc, j-2, i  );
2987   b13 = tex2D(tex_img_inc, j+3, i  );
2988   //minmax aux extremites
2989   minmax5(&a9,&a10,&a11,&a12,&a13);
2990   minmax5(&b9,&b10,&b11,&b12,&b13);
2991
2992   //chargement valeur suivante (24)
2993   a13 = tex2D(tex_img_inc, j-2, i+1);
2994   b13 = tex2D(tex_img_inc, j+3, i+1);
2995   //minmax aux extremites
2996   minmax4(&a10,&a11,&a12,&a13);
2997   minmax4(&b10,&b11,&b12,&b13);
2998
2999   //chargement valeur suivante (25)
3000   a13 = tex2D(tex_img_inc, j-2, i+2);
3001   b13 = tex2D(tex_img_inc, j+3, i+2);
3002   //minmax aux extremites
3003   minmax3(&a11,&a12,&a13);
3004   minmax3(&b11,&b12,&b13);
3005
3006   //median au milieu !
3007   output[ __mul24(i, j_dim) +j ] = a12 ;
3008   output[ __mul24(i, j_dim) +j+1 ] = b12 ;
3009
3010 }
3011
3012
3013 __global__ void kernel_median7( unsigned short*output, int i_dim, int j_dim)
3014 {
3015
3016   // coordonnees absolues du point
3017   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
3018   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
3019
3020   /**************************************************************************
3021    *             tri(s)
3022    **************************************************************************/
3023   int a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,a25 ;
3024
3025   /****************************************
3026    * les 26 premieres valeurs
3027    ****************************************/
3028   //premiere ligne
3029   a0  = tex2D(tex_img_ins, j-2, i-3) ;
3030   a1  = tex2D(tex_img_ins, j-1, i-3) ;
3031   a2  = tex2D(tex_img_ins, j  , i-3) ;
3032   a3  = tex2D(tex_img_ins, j+1, i-3) ;
3033   a4  = tex2D(tex_img_ins, j+2, i-3) ;
3034   a5  = tex2D(tex_img_ins, j+3, i-3) ;
3035   //deuxieme ligne
3036   a6  = tex2D(tex_img_ins, j-2, i-2) ;
3037   a7  = tex2D(tex_img_ins, j-1, i-2) ;
3038   a8  = tex2D(tex_img_ins, j  , i-2) ;
3039   a9  = tex2D(tex_img_ins, j+1, i-2) ;
3040   a10 = tex2D(tex_img_ins, j+2, i-2) ;
3041   a11 = tex2D(tex_img_ins, j+3, i-2) ;
3042   //troisieme ligne
3043   a12 = tex2D(tex_img_ins, j-2, i-1) ;
3044   a13 = tex2D(tex_img_ins, j-1, i-1) ;
3045   a14 = tex2D(tex_img_ins, j  , i-1) ;
3046   a15 = tex2D(tex_img_ins, j+1, i-1) ;
3047   a16 = tex2D(tex_img_ins, j+2, i-1) ;
3048   a17 = tex2D(tex_img_ins, j+3, i-1) ;
3049   //quatrieme ligne
3050   a18 = tex2D(tex_img_ins, j-2, i  ) ;
3051   a19 = tex2D(tex_img_ins, j-1, i  ) ;
3052   a20 = tex2D(tex_img_ins, j  , i  ) ;
3053   a21 = tex2D(tex_img_ins, j+1, i  ) ;
3054   a22 = tex2D(tex_img_ins, j+2, i  ) ;
3055   a23 = tex2D(tex_img_ins, j+3, i  ) ;
3056   //cinquieme ligne
3057   a24 = tex2D(tex_img_ins, j-2, i+1) ;
3058   a25 = tex2D(tex_img_ins, j-1, i+1) ;
3059
3060   //min max aux extremites
3061   minmax26(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3062
3063   //chargement valeurs suivante (26)
3064   a25 = tex2D(tex_img_ins, j , i+1);
3065   //minmax aux extremites
3066   minmax25(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3067
3068   //chargement valeur suivante (27)
3069   a25 = tex2D(tex_img_ins, j+1, i+1);
3070   //minmax aux extremites
3071   minmax24(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3072
3073   //chargement valeur suivante (28)
3074   a25 = tex2D(tex_img_ins, j+2, i+1);
3075   //minmax aux extremites
3076   minmax23(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3077
3078   //chargement valeur suivante (29)
3079   a25 = tex2D(tex_img_ins, j+3, i+1);
3080   //minmax aux extremites
3081   minmax22(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3082
3083   //chargement valeur suivante (30)
3084   a25 = tex2D(tex_img_ins, j-2, i+2);
3085   //minmax aux extremites
3086   minmax21(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3087
3088   //chargement valeur suivante (31)
3089   a25 = tex2D(tex_img_ins, j-1, i+2);
3090   //minmax aux extmites
3091   minmax20(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3092
3093   //chargement valeur suivante (32)
3094   a25 = tex2D(tex_img_ins, j  , i+2);
3095   //minmax aux extmites
3096   minmax19(&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3097
3098   //chargement valeur suivante (33)
3099   a25 = tex2D(tex_img_ins, j+1, i+2);
3100   //minmax aux extmites
3101   minmax18(&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3102
3103   //chargement valeur suivante (34)
3104   a25 = tex2D(tex_img_ins, j+2, i+2);
3105   //minmax aux extmites
3106   minmax17(&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3107
3108   //chargement valeur suivante (35)
3109   a25 = tex2D(tex_img_ins, j+3, i+2);
3110   //minmax aux extmites
3111   minmax16(&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3112
3113   //chargement valeur suivante (36)
3114   a25 = tex2D(tex_img_ins, j-2, i+3);
3115   //minmax aux extmites
3116   minmax15(&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3117
3118   //chargement valeur suivante (37)
3119   a25 = tex2D(tex_img_ins, j-1, i+3);
3120   //minmax aux extmites
3121
3122   minmax14(&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3123   //chargement valeur suivante (38)
3124   a25 = tex2D(tex_img_ins, j  , i+3);
3125   //minmax aux extmites
3126   minmax13(&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3127
3128   //chargement valeur suivante (39)
3129   a25 = tex2D(tex_img_ins, j+1, i+3);
3130   //minmax aux extmites
3131   minmax12(&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3132
3133   //chargement valeur suivante (40)
3134   a25 = tex2D(tex_img_ins, j+2, i+3);
3135   //minmax aux extmites
3136   minmax11(&a15, &a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3137
3138   //chargement valeur suivante (41)
3139   a25 = tex2D(tex_img_ins, j+3, i+3);
3140   //minmax aux extmites
3141   minmax10(&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3142
3143   //chargement valeurs suivantes
3144   a25 = tex2D(tex_img_ins, j-3, i+3);
3145   //minmax aux extremites
3146   minmax9(&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3147
3148   //chargement valeurs suivantes
3149   a25 = tex2D(tex_img_ins, j-3, i+2);
3150   //minmax aux extremites
3151   minmax8(&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3152
3153   //chargement valeurs suivantes
3154   a25 = tex2D(tex_img_ins, j-3, i+1);
3155   //minmax aux extremites
3156   minmax7(&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3157
3158   //chargement valeurs suivantes
3159   a25 = tex2D(tex_img_ins, j-3, i);
3160   //minmax aux extremites
3161   minmax6(&a20,&a21,&a22,&a23,&a24,&a25);
3162
3163   //chargement valeurs suivantes
3164   a25 = tex2D(tex_img_ins, j-3, i-1);
3165   //minmax aux extremites
3166   minmax5(&a21,&a22,&a23,&a24,&a25);
3167
3168   //chargement valeurs suivantes
3169   a25 = tex2D(tex_img_ins, j-3, i-2);
3170   //minmax aux extremites
3171   minmax4(&a22,&a23,&a24,&a25);
3172
3173   //chargement valeurs suivantes
3174   a25 = tex2D(tex_img_ins, j-3, i-3);
3175   //minmax aux extremites
3176   minmax3(&a23,&a24,&a25);
3177
3178   //medians au milieu !
3179   output[ __mul24(i, j_dim) +j ] = a24 ;
3180
3181 }
3182
3183
3184
3185 __global__ void kernel_median7_2pix( unsigned short*output, int i_dim, int j_dim)
3186 {
3187
3188   // coordonnees absolues du point
3189   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
3190   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
3191
3192   /**************************************************************************
3193    *             tri(s)
3194    **************************************************************************/
3195   int a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,a25 ;
3196   int b17,b18,b19,b20,b21,b22,b23,b24,b25;
3197
3198   /****************************************
3199    * les 26 premieres valeurs
3200    ****************************************/
3201   //premiere ligne
3202   a0  = tex2D(tex_img_ins, j-2, i-3) ;
3203   a1  = tex2D(tex_img_ins, j-1, i-3) ;
3204   a2  = tex2D(tex_img_ins, j  , i-3) ;
3205   a3  = tex2D(tex_img_ins, j+1, i-3) ;
3206   a4  = tex2D(tex_img_ins, j+2, i-3) ;
3207   a5  = tex2D(tex_img_ins, j+3, i-3) ;
3208   //deuxieme ligne
3209   a6  = tex2D(tex_img_ins, j-2, i-2) ;
3210   a7  = tex2D(tex_img_ins, j-1, i-2) ;
3211   a8  = tex2D(tex_img_ins, j  , i-2) ;
3212   a9  = tex2D(tex_img_ins, j+1, i-2) ;
3213   a10 = tex2D(tex_img_ins, j+2, i-2) ;
3214   a11 = tex2D(tex_img_ins, j+3, i-2) ;
3215   //troisieme ligne
3216   a12 = tex2D(tex_img_ins, j-2, i-1) ;
3217   a13 = tex2D(tex_img_ins, j-1, i-1) ;
3218   a14 = tex2D(tex_img_ins, j  , i-1) ;
3219   a15 = tex2D(tex_img_ins, j+1, i-1) ;
3220   a16 = tex2D(tex_img_ins, j+2, i-1) ;
3221   a17 = tex2D(tex_img_ins, j+3, i-1) ;
3222   //quatrieme ligne
3223   a18 = tex2D(tex_img_ins, j-2, i  ) ;
3224   a19 = tex2D(tex_img_ins, j-1, i  ) ;
3225   a20 = tex2D(tex_img_ins, j  , i  ) ;
3226   a21 = tex2D(tex_img_ins, j+1, i  ) ;
3227   a22 = tex2D(tex_img_ins, j+2, i  ) ;
3228   a23 = tex2D(tex_img_ins, j+3, i  ) ;
3229   //cinquieme ligne
3230   a24 = tex2D(tex_img_ins, j-2, i+1) ;
3231   a25 = tex2D(tex_img_ins, j-1, i+1) ;
3232
3233   //min max aux extremites
3234   minmax26(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3235
3236   //chargement valeurs suivante (26)
3237   a25 = tex2D(tex_img_ins, j , i+1);
3238   //minmax aux extremites
3239   minmax25(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3240
3241   //chargement valeur suivante (27)
3242   a25 = tex2D(tex_img_ins, j+1, i+1);
3243   //minmax aux extremites
3244   minmax24(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3245
3246   //chargement valeur suivante (28)
3247   a25 = tex2D(tex_img_ins, j+2, i+1);
3248   //minmax aux extremites
3249   minmax23(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3250
3251   //chargement valeur suivante (29)
3252   a25 = tex2D(tex_img_ins, j+3, i+1);
3253   //minmax aux extremites
3254   minmax22(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3255
3256   //chargement valeur suivante (30)
3257   a25 = tex2D(tex_img_ins, j-2, i+2);
3258   //minmax aux extremites
3259   minmax21(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3260
3261   //chargement valeur suivante (31)
3262   a25 = tex2D(tex_img_ins, j-1, i+2);
3263   //minmax aux extmites
3264   minmax20(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3265
3266   //chargement valeur suivante (32)
3267   a25 = tex2D(tex_img_ins, j  , i+2);
3268   //minmax aux extmites
3269   minmax19(&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3270
3271   //chargement valeur suivante (33)
3272   a25 = tex2D(tex_img_ins, j+1, i+2);
3273   //minmax aux extmites
3274   minmax18(&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3275
3276   //chargement valeur suivante (34)
3277   a25 = tex2D(tex_img_ins, j+2, i+2);
3278   //minmax aux extmites
3279   minmax17(&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3280
3281   //chargement valeur suivante (35)
3282   a25 = tex2D(tex_img_ins, j+3, i+2);
3283   //minmax aux extmites
3284   minmax16(&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3285
3286   //chargement valeur suivante (36)
3287   a25 = tex2D(tex_img_ins, j-2, i+3);
3288   //minmax aux extmites
3289   minmax15(&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3290
3291   //chargement valeur suivante (37)
3292   a25 = tex2D(tex_img_ins, j-1, i+3);
3293   //minmax aux extmites
3294
3295   minmax14(&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3296   //chargement valeur suivante (38)
3297   a25 = tex2D(tex_img_ins, j  , i+3);
3298   //minmax aux extmites
3299   minmax13(&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3300
3301   //chargement valeur suivante (39)
3302   a25 = tex2D(tex_img_ins, j+1, i+3);
3303   //minmax aux extmites
3304   minmax12(&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3305
3306   //chargement valeur suivante (40)
3307   a25 = tex2D(tex_img_ins, j+2, i+3);
3308   //minmax aux extmites
3309   minmax11(&a15, &a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3310
3311   //chargement valeur suivante (41)
3312   a25 = tex2D(tex_img_ins, j+3, i+3);
3313   //minmax aux extmites
3314   minmax10(&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3315
3316   // fin des pixels voisins communs deux pixels centraux
3317   b17=a17; b18=a18; b19=a19; b20=a20; b21=a21; b22=a22; b23=a23; b24=a24; b25=a25;
3318
3319   //chargement valeurs suivantes
3320   a25 = tex2D(tex_img_ins, j-3, i+3);
3321   b25 = tex2D(tex_img_ins, j+4, i+3);
3322   //minmax aux extremites
3323   minmax9(&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3324   minmax9(&b17,&b18,&b19,&b20,&b21,&b22,&b23,&b24,&b25);
3325
3326   //chargement valeurs suivantes
3327   a25 = tex2D(tex_img_ins, j-3, i+2);
3328   b25 = tex2D(tex_img_ins, j+4, i+2);
3329   //minmax aux extremites
3330   minmax8(&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3331   minmax8(&b18,&b19,&b20,&b21,&b22,&b23,&b24,&b25);
3332
3333   //chargement valeurs suivantes
3334   a25 = tex2D(tex_img_ins, j-3, i+1);
3335   b25 = tex2D(tex_img_ins, j+4, i+1);
3336   //minmax aux extremites
3337   minmax7(&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3338   minmax7(&b19,&b20,&b21,&b22,&b23,&b24,&b25);
3339
3340   //chargement valeurs suivantes
3341   a25 = tex2D(tex_img_ins, j-3, i);
3342   b25 = tex2D(tex_img_ins, j+4, i);
3343   //minmax aux extremites
3344   minmax6(&a20,&a21,&a22,&a23,&a24,&a25);
3345   minmax6(&b20,&b21,&b22,&b23,&b24,&b25);
3346
3347   //chargement valeurs suivantes
3348   a25 = tex2D(tex_img_ins, j-3, i-1);
3349   b25 = tex2D(tex_img_ins, j+4, i-1);
3350   //minmax aux extremites
3351   minmax5(&a21,&a22,&a23,&a24,&a25);
3352   minmax5(&b21,&b22,&b23,&b24,&b25);
3353
3354   //chargement valeurs suivantes
3355   a25 = tex2D(tex_img_ins, j-3, i-2);
3356   b25 = tex2D(tex_img_ins, j+4, i-2);
3357   //minmax aux extremites
3358   minmax4(&a22,&a23,&a24,&a25);
3359   minmax4(&b22,&b23,&b24,&b25);
3360
3361   //chargement valeurs suivantes
3362   a25 = tex2D(tex_img_ins, j-3, i-3);
3363   b25 = tex2D(tex_img_ins, j+4, i-3);
3364   //minmax aux extremites
3365   minmax3(&a23,&a24,&a25);
3366   minmax3(&b23,&b24,&b25);
3367
3368   //medians au milieu !
3369   output[ __mul24(i, j_dim) +j ] = a24 ;
3370   output[ __mul24(i, j_dim) +j+1 ] = b24 ;
3371
3372 }
3373
3374
3375 __global__ void kernel_median7_2pix( unsigned char  *output, int i_dim, int j_dim)
3376 {
3377
3378   // coordonnees absolues du point
3379   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
3380   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
3381
3382   /**************************************************************************
3383    *             tri(s)
3384    **************************************************************************/
3385   int a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,a25 ;
3386   int b17,b18,b19,b20,b21,b22,b23,b24,b25;
3387
3388   /****************************************
3389    * les 26 premieres valeurs
3390    ****************************************/
3391   //premiere ligne
3392   a0  = tex2D(tex_img_inc, j-2, i-3) ;
3393   a1  = tex2D(tex_img_inc, j-1, i-3) ;
3394   a2  = tex2D(tex_img_inc, j  , i-3) ;
3395   a3  = tex2D(tex_img_inc, j+1, i-3) ;
3396   a4  = tex2D(tex_img_inc, j+2, i-3) ;
3397   a5  = tex2D(tex_img_inc, j+3, i-3) ;
3398   //deuxieme ligne
3399   a6  = tex2D(tex_img_inc, j-2, i-2) ;
3400   a7  = tex2D(tex_img_inc, j-1, i-2) ;
3401   a8  = tex2D(tex_img_inc, j  , i-2) ;
3402   a9  = tex2D(tex_img_inc, j+1, i-2) ;
3403   a10 = tex2D(tex_img_inc, j+2, i-2) ;
3404   a11 = tex2D(tex_img_inc, j+3, i-2) ;
3405   //troisieme ligne
3406   a12 = tex2D(tex_img_inc, j-2, i-1) ;
3407   a13 = tex2D(tex_img_inc, j-1, i-1) ;
3408   a14 = tex2D(tex_img_inc, j  , i-1) ;
3409   a15 = tex2D(tex_img_inc, j+1, i-1) ;
3410   a16 = tex2D(tex_img_inc, j+2, i-1) ;
3411   a17 = tex2D(tex_img_inc, j+3, i-1) ;
3412   //quatrieme ligne
3413   a18 = tex2D(tex_img_inc, j-2, i  ) ;
3414   a19 = tex2D(tex_img_inc, j-1, i  ) ;
3415   a20 = tex2D(tex_img_inc, j  , i  ) ;
3416   a21 = tex2D(tex_img_inc, j+1, i  ) ;
3417   a22 = tex2D(tex_img_inc, j+2, i  ) ;
3418   a23 = tex2D(tex_img_inc, j+3, i  ) ;
3419   //cinquieme ligne
3420   a24 = tex2D(tex_img_inc, j-2, i+1) ;
3421   a25 = tex2D(tex_img_inc, j-1, i+1) ;
3422
3423   //min max aux extremites
3424   minmax26(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3425
3426   //chargement valeurs suivante (26)
3427   a25 = tex2D(tex_img_inc, j , i+1);
3428   //minmax aux extremites
3429   minmax25(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3430
3431   //chargement valeur suivante (27)
3432   a25 = tex2D(tex_img_inc, j+1, i+1);
3433   //minmax aux extremites
3434   minmax24(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3435
3436   //chargement valeur suivante (28)
3437   a25 = tex2D(tex_img_inc, j+2, i+1);
3438   //minmax aux extremites
3439   minmax23(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3440
3441   //chargement valeur suivante (29)
3442   a25 = tex2D(tex_img_inc, j+3, i+1);
3443   //minmax aux extremites
3444   minmax22(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3445
3446   //chargement valeur suivante (30)
3447   a25 = tex2D(tex_img_inc, j-2, i+2);
3448   //minmax aux extremites
3449   minmax21(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3450
3451   //chargement valeur suivante (31)
3452   a25 = tex2D(tex_img_inc, j-1, i+2);
3453   //minmax aux extmites
3454   minmax20(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3455
3456   //chargement valeur suivante (32)
3457   a25 = tex2D(tex_img_inc, j  , i+2);
3458   //minmax aux extmites
3459   minmax19(&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3460
3461   //chargement valeur suivante (33)
3462   a25 = tex2D(tex_img_inc, j+1, i+2);
3463   //minmax aux extmites
3464   minmax18(&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3465
3466   //chargement valeur suivante (34)
3467   a25 = tex2D(tex_img_inc, j+2, i+2);
3468   //minmax aux extmites
3469   minmax17(&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3470
3471   //chargement valeur suivante (35)
3472   a25 = tex2D(tex_img_inc, j+3, i+2);
3473   //minmax aux extmites
3474   minmax16(&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3475
3476   //chargement valeur suivante (36)
3477   a25 = tex2D(tex_img_inc, j-2, i+3);
3478   //minmax aux extmites
3479   minmax15(&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3480
3481   //chargement valeur suivante (37)
3482   a25 = tex2D(tex_img_inc, j-1, i+3);
3483   //minmax aux extmites
3484
3485   minmax14(&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3486   //chargement valeur suivante (38)
3487   a25 = tex2D(tex_img_inc, j  , i+3);
3488   //minmax aux extmites
3489   minmax13(&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3490
3491   //chargement valeur suivante (39)
3492   a25 = tex2D(tex_img_inc, j+1, i+3);
3493   //minmax aux extmites
3494   minmax12(&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3495
3496   //chargement valeur suivante (40)
3497   a25 = tex2D(tex_img_inc, j+2, i+3);
3498   //minmax aux extmites
3499   minmax11(&a15, &a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3500
3501   //chargement valeur suivante (41)
3502   a25 = tex2D(tex_img_inc, j+3, i+3);
3503   //minmax aux extmites
3504   minmax10(&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3505
3506   // fin des pixels voisins communs deux pixels centraux
3507   b17=a17; b18=a18; b19=a19; b20=a20; b21=a21; b22=a22; b23=a23; b24=a24; b25=a25;
3508
3509   //chargement valeurs suivantes
3510   a25 = tex2D(tex_img_inc, j-3, i+3);
3511   b25 = tex2D(tex_img_inc, j+4, i+3);
3512   //minmax aux extremites
3513   minmax9(&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3514   minmax9(&b17,&b18,&b19,&b20,&b21,&b22,&b23,&b24,&b25);
3515
3516   //chargement valeurs suivantes
3517   a25 = tex2D(tex_img_inc, j-3, i+2);
3518   b25 = tex2D(tex_img_inc, j+4, i+2);
3519   //minmax aux extremites
3520   minmax8(&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3521   minmax8(&b18,&b19,&b20,&b21,&b22,&b23,&b24,&b25);
3522
3523   //chargement valeurs suivantes
3524   a25 = tex2D(tex_img_inc, j-3, i+1);
3525   b25 = tex2D(tex_img_inc, j+4, i+1);
3526   //minmax aux extremites
3527   minmax7(&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3528   minmax7(&b19,&b20,&b21,&b22,&b23,&b24,&b25);
3529
3530   //chargement valeurs suivantes
3531   a25 = tex2D(tex_img_inc, j-3, i);
3532   b25 = tex2D(tex_img_inc, j+4, i);
3533   //minmax aux extremites
3534   minmax6(&a20,&a21,&a22,&a23,&a24,&a25);
3535   minmax6(&b20,&b21,&b22,&b23,&b24,&b25);
3536
3537   //chargement valeurs suivantes
3538   a25 = tex2D(tex_img_inc, j-3, i-1);
3539   b25 = tex2D(tex_img_inc, j+4, i-1);
3540   //minmax aux extremites
3541   minmax5(&a21,&a22,&a23,&a24,&a25);
3542   minmax5(&b21,&b22,&b23,&b24,&b25);
3543
3544   //chargement valeurs suivantes
3545   a25 = tex2D(tex_img_inc, j-3, i-2);
3546   b25 = tex2D(tex_img_inc, j+4, i-2);
3547   //minmax aux extremites
3548   minmax4(&a22,&a23,&a24,&a25);
3549   minmax4(&b22,&b23,&b24,&b25);
3550
3551   //chargement valeurs suivantes
3552   a25 = tex2D(tex_img_inc, j-3, i-3);
3553   b25 = tex2D(tex_img_inc, j+4, i-3);
3554   //minmax aux extremites
3555   minmax3(&a23,&a24,&a25);
3556   minmax3(&b23,&b24,&b25);
3557
3558   //medians au milieu !
3559   output[ __mul24(i, j_dim) +j ] = a24 ;
3560   output[ __mul24(i, j_dim) +j+1 ] = b24 ;
3561
3562 }
3563
3564
3565
3566 /*****************************************************************************
3567  *      median generic shared mem + forgetfull
3568  *****************************************************************************/
3569 __global__ void kernel_medianForgetRSH( unsigned short * output, int i_dim, int j_dim, int r)
3570 {
3571   // coordonnees absolues du point
3572   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
3573   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
3574   unsigned short ic, jc ;
3575   unsigned short Nreg = ((2*r+1)*(2*r+1))/2 + 2 ;
3576
3577   int bdimX = blockDim.x; //<<1 ; // pour le cas à 2 pixels par thread
3578   int tidX = threadIdx.x; //<<1 ;
3579
3580   // chargement en smem
3581   int idrow = threadIdx.y*(blockDim.x+2*r) ;
3582
3583   extern __shared__ int medRoi[];
3584
3585   // bloc 0 (en haut à gauche)
3586   medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-r, i-r) ;
3587   // bloc 1 (en haut à droite)...
3588   if ( threadIdx.x < 2*r ) //...ou plutot ce qu'il en manque
3589         medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-r, i-r ) ;
3590   // bloc 2 ( en bas à gauche)
3591   if ( threadIdx.y < 2*r )
3592         {
3593           idrow = (threadIdx.y+blockDim.y)*(blockDim.x+2*r) ;
3594           medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-r, i+blockDim.y-r ) ;
3595           //bloc 4 ( en bas à doite )...
3596           if ( threadIdx.x < 2*r ) //...ou ce qu'il en manque
3597                 medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-r, i+blockDim.y-r ) ;
3598         }
3599   __syncthreads() ;
3600
3601   // remplissage du vecteur de tri minmax
3602   unsigned short vect[8066] ;
3603   int Freg=Nreg ;
3604   for (ic=0; ic<2*r+1; ic++)
3605         {
3606           for (jc=0; jc<2*r+1; jc++)
3607                 {
3608                   if ( ic*(2*r+1)+jc < Nreg )
3609                         {
3610                           vect[ ic*(2*r+1)+jc ] = medRoi[ (threadIdx.y+ic)*(bdimX+2*r)+ (tidX+jc) ] ;
3611                         } else
3612                         {
3613                           minmaxN(vect, Freg--) ;
3614                           vect[ Nreg-1 ] = medRoi[ (threadIdx.y+ic)*(bdimX+2*r)+ (tidX+jc) ] ;
3615                         }
3616                 }
3617         }
3618   minmax3(&vect[Nreg-3], &vect[Nreg-2], &vect[Nreg-1])
3619
3620   //medRoi[ (threadIdx.y+ic)*(bdimX+L-1)+ (tidX+jc) ]
3621
3622         output[ __mul24(i, j_dim) +j ] = vect[ Nreg-2 ];
3623 }
3624
3625
3626 /*****************************************************************************
3627  *      median generic shared mem + forgetfull
3628  *****************************************************************************/
3629 __global__ void kernel_medianForgetR( unsigned short * output, int i_dim, int j_dim, int r)
3630 {
3631   // coordonnees absolues du point
3632   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
3633   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
3634   unsigned short ic, jc ;
3635   unsigned short Nreg = ((2*r+1)*(2*r+1))/2 + 2 ;
3636
3637   // remplissage du vecteur de tri minmax
3638   unsigned short vect[8066] ;
3639   int Freg=Nreg ;
3640   for (ic=0; ic<2*r+1; ic++)
3641         {
3642           for (jc=0; jc<2*r+1; jc++)
3643                 {
3644                   if ( ic*(2*r+1)+jc < Nreg )
3645                         {
3646                           vect[ ic*(2*r+1)+jc ] = tex2D(tex_img_ins, j-r+jc, i-r+ic) ;
3647                         } else
3648                         {
3649                           minmaxN(vect, Freg--) ;
3650                           vect[ Nreg-1 ] = tex2D(tex_img_ins, j-r+jc, i-r+ic) ;
3651                         }
3652                 }
3653         }
3654   minmax3(&vect[Nreg-3], &vect[Nreg-2], &vect[Nreg-1])
3655
3656   //medRoi[ (threadIdx.y+ic)*(bdimX+L-1)+ (tidX+jc) ]
3657
3658         output[ __mul24(i, j_dim) +j ] = vect[ Nreg-2 ];
3659 }
3660
3661
3662
3663 /**********************************************************************************
3664  *    MEDIAN PSEUDO-SEPARABLE POUR GRANDS NOYAUX
3665  **********************************************************************************/
3666 __global__ void kernel_medianSepR( unsigned short *output, int i_dim, int j_dim, int r)
3667 {
3668
3669   int idc, val, min, max, inf, egal, sup, mxinf, minsup, estim ;
3670
3671   //coordonnées ds le bloc
3672   int ib = threadIdx.y ;
3673   int jb = threadIdx.x ;
3674   //int idx_h = __mul24(ib+r,blockDim.x) + jb ;   // index pixel deans shmem (bloc+halo)
3675   //int offset = __mul24(blockDim.x,r) ;
3676
3677   // coordonnees absolues du point
3678   int j = __mul24(blockIdx.x,blockDim.x) + jb ;
3679   int i = __mul24(blockIdx.y,blockDim.y) + ib ;
3680
3681   extern __shared__ int buff[] ;
3682   /***********************************************************************************
3683    *              CHARGEMENT DATA EN SHARED MEM
3684    ***********************************************************************************/
3685   for (idc = 0 ; idc < (2*(blockDim.y+r)-1)/blockDim.y; idc++)
3686         {
3687           if (idc*blockDim.y +ib < i_dim)
3688                 buff[ (idc*blockDim.y +ib)*blockDim.x + jb ]    = tex2D(tex_img_ins, j, i-r+ idc*blockDim.y) ;
3689         }
3690
3691   __syncthreads() ;
3692   /**********************************************************************************************
3693    *               TRI VERTICAL par algo TORBEN MOGENSEN
3694    *          (a little bit slow but saves memory => faster !)
3695    **********************************************************************************************/
3696   min = max = buff[ ib*blockDim.x +jb] ;
3697
3698   for (idc= 0 ; idc< 2*r+1 ; idc++ )
3699         {
3700           val = buff[ __mul24(ib+idc, blockDim.x) +jb ] ;
3701           if ( val  < min ) min = val ;
3702           if ( val  > max ) max = val ;
3703         }
3704
3705   while (1)
3706         {
3707           estim = (min+max)/2 ;
3708           inf = sup = egal = 0  ;
3709           mxinf = min ;
3710           minsup= max ;
3711           for (idc =0; idc< 2*r+1 ; idc++)
3712                 {
3713                   val = buff[ __mul24(ib+idc, blockDim.x) +jb ] ;
3714                   if( val < estim )
3715                         {
3716                           inf++;
3717                           if( val > mxinf) mxinf = val ;
3718                         } else if (val > estim)
3719                         {
3720                           sup++;
3721                           if( val < minsup) minsup = val ;
3722                         } else egal++ ;
3723                 }
3724           if ( (inf <= (r+1))&&(sup <=(r+1)) ) break ;
3725           else if (inf>sup) max = mxinf ;
3726           else min = minsup ;
3727           }
3728
3729   if ( inf >= r+1 ) val = mxinf ;
3730   else if (inf+egal >= r+1) val = estim ;
3731   else val = minsup ;
3732
3733   output[ __mul24(j, i_dim)  +i  ] =  val ;
3734 }
3735
3736
3737 /**
3738  *
3739  * correction de staircase
3740  */
3741 __global__ void kernel_staircase_reduc3(unsigned int * img_out, unsigned int L, unsigned int H)
3742 {
3743   // coordonnees du point dans le bloc
3744   //unsigned int iib = threadIdx.x ;
3745   //unsigned int jib = threadIdx.y ;
3746   // coordonnees du point dans l'image
3747   unsigned int y = blockIdx.y*blockDim.y + threadIdx.y;
3748   unsigned int x = blockIdx.x*blockDim.x + threadIdx.x;
3749
3750   int a, b, c, d, e, f, g, h, i ;            // gl des voisins
3751   float wa, wb, wc, wd, we, wf, wg, wh, wi ; // poids
3752   float S1, S2, S11, S22, S12, S0, Sx, Sx1, Sx2 ;
3753   float c1,c2 ;
3754
3755   // chargement des valeurs GL des pixels du voisinage
3756   a = tex2D(tex_img_estim, x-1, y-1) ;
3757   b = tex2D(tex_img_estim, x  , y-1) ;
3758   c = tex2D(tex_img_estim, x+1, y-1) ;
3759   d = tex2D(tex_img_estim, x-1, y  ) ;
3760   e = tex2D(tex_img_estim, x  , y  ) ;
3761   f = tex2D(tex_img_estim, x+1, y  ) ;
3762   g = tex2D(tex_img_estim, x-1, y+1) ;
3763   h = tex2D(tex_img_estim, x  , y+1) ;
3764   i = tex2D(tex_img_estim, x+1, y+1) ;
3765
3766   wa = tex1D(tex_noyau, abs(a-e)) ;
3767   wb = tex1D(tex_noyau, abs(b-e)) ;
3768   wc = tex1D(tex_noyau, abs(c-e)) ;
3769   wd = tex1D(tex_noyau, abs(d-e)) ;
3770   we = 1 ;
3771   wf = tex1D(tex_noyau, abs(f-e)) ;
3772   wg = tex1D(tex_noyau, abs(g-e)) ;
3773   wh = tex1D(tex_noyau, abs(h-e)) ;
3774   wi = tex1D(tex_noyau, abs(i-e)) ;
3775
3776
3777   //construction des elements du systeme lineaire
3778   S0  = wa+wb+wc+wd+we+wf+wg+wh+wi ;
3779   S1  = wc+wf+wi-wa-wd-wg ;
3780   S2  = wg+wh+wi-wa-wb-wc ;
3781   Sx  = wa*a + wb*b + wc*c + wd*d + we*e + wf*f + wg*g + wh*h + wi*i ;
3782   Sx1 = wc*c + wf*f + wi*i - wa*a - wd*d - wg*g ;
3783   Sx2 = wg*g + wh*h + wi*i - wa*a - wb*b - wc*c ;
3784   S11 = wc+wf+wi+wa+wd+wg ;
3785   S22 = wg+wh+wi+wa+wb+wc ;
3786   S12 = wa + wi -wc -wg ;
3787
3788   c1 = S22*(S11*Sx-S1*Sx1) + S12*(S1*Sx2-S12*Sx) + S2*(S12*Sx1-S11*Sx2) ;
3789   c2 = S22*(S11*S0-S1*S1) + S12*(S2*S1-S12*S0) + S2*(S1*S12-S2*S11) ;
3790   img_out[y*L + x] = c1/c2 ;
3791 }
3792
3793