THESE/Chapters/chapter6/code/convoSepoptimH.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, base ;
 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           base = __umul24(ic,L) ;
 459           pix = tex2D(tex_img_inc, j+1, i-1+ic) ;
 460           outval0 += masque[ base +2 ]*pix ;
 461           outval1 += masque[ base +1 ]*pix ;
 462           outval2 += masque[ base    ]*pix ;
 463           pix = tex2D(tex_img_inc, j+2, i-1+ic) ;
 464           outval1 += masque[ base +2 ]*pix ;
 465           outval2 += masque[ base +1 ]*pix ;
 466           outval3 += masque[ base    ]*pix ;
 467           pix = tex2D(tex_img_inc, j+3, i-1+ic) ;
 468           outval2 += masque[ base +2 ]*pix ;
 469           outval3 += masque[ base +1 ]*pix ;
 470           outval4 += masque[ base    ]*pix ;
 471           pix = tex2D(tex_img_inc, j+4, i-1+ic) ;
 472           outval3 += masque[ base +2 ]*pix ;
 473           outval4 += masque[ base +1 ]*pix ;
 474           outval5 += masque[ base    ]*pix ;
 475           pix = tex2D(tex_img_inc, j+5, i-1+ic) ;
 476           outval4 += masque[ base +2 ]*pix ;
 477           outval5 += masque[ base +1 ]*pix ;
 478           outval6 += masque[ base    ]*pix ;
 479           pix = tex2D(tex_img_inc, j+6, i-1+ic) ;
 480           outval5 += masque[ base +2 ]*pix ;
 481           outval6 += masque[ base +1 ]*pix ;
 482           outval7 += masque[ base    ]*pix ;
 483
 484           pix = tex2D(tex_img_inc, j, i-1+ic) ;
 485           outval0 += masque[ base +1 ]*pix ;
 486           outval1 += masque[ base    ]*pix ;
 487           pix = tex2D(tex_img_inc, j-1, i-1+ic) ;
 488           outval0 += masque[ base  ]*pix ;
 489
 490           pix = tex2D(tex_img_inc, j+7, i-1+ic) ;
 491           outval6 += masque[ base +2 ]*pix ;
 492           outval7 += masque[ base +1 ]*pix ;
 493           pix = tex2D(tex_img_inc, j+8, i-1+ic) ;
 494           outval7 += masque[ base +2 ]*pix ;
 495         }
 496
 497   // les 2 pixels hauts
 498   base = __umul24(i, j_dim) + j ;
 499   output[ base++ ] = outval0 ;
 500   output[ base++ ] = outval1 ;
 501   output[ base++ ] = outval2 ;
 502   output[ base++ ] = outval3;
 503   output[ base++ ] = outval4;
 504   output[ base++ ] = outval5;
 505   output[ base++ ] = outval6 ;
 506   output[ base++ ] = outval7 ;
 507
 508 }
 509
 510
 511
 512
 513
 514 // convolution non separable 8 pixels/thread
 515 // image en texture et noyau en mem constante
 516 __global__ void kernel_convoGene8x8pL5( unsigned char  *output, int j_dim )
 517 {
 518   int ic, base ;
 519   int L=5 ;
 520   unsigned char pix ;
 521   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
 522   float outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
 523
 524   // coordonnees absolues du point de base en haut a gauche
 525   int j = ( __umul24( blockIdx.x, blockDim.x) + threadIdx.x)<< 3 ;
 526   int i = ( __umul24( blockIdx.y, blockDim.y) + threadIdx.y) ;
 527
 528   //zone commune
 529   //forme generique contrib pixel => masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic)
 530   for (ic=0 ; ic<L ; ic++)
 531         {
 532           base = __umul24(ic,L) ;
 533           pix = tex2D(tex_img_inc, j+2, i-2+ic) ;
 534           outval0 += masque[ base +4 ]*pix ;
 535           outval1 += masque[ base +3 ]*pix ;
 536           outval2 += masque[ base +2 ]*pix ;
 537           outval3 += masque[ base +1 ]*pix ;
 538           outval4 += masque[ base    ]*pix ;
 539           pix = tex2D(tex_img_inc, j+3, i-2+ic) ;
 540           outval1 += masque[ base +4 ]*pix ;
 541           outval2 += masque[ base +3 ]*pix ;
 542           outval3 += masque[ base +2 ]*pix ;
 543           outval4 += masque[ base +1 ]*pix ;
 544           outval5 += masque[ base    ]*pix ;
 545           pix = tex2D(tex_img_inc, j+4, i-2+ic) ;
 546           outval2 += masque[ base +4 ]*pix ;
 547           outval3 += masque[ base +3 ]*pix ;
 548           outval4 += masque[ base +2 ]*pix ;
 549           outval5 += masque[ base +1 ]*pix ;
 550           outval6 += masque[ base    ]*pix ;
 551           pix = tex2D(tex_img_inc, j+5, i-2+ic) ;
 552           outval3 += masque[ base +4 ]*pix ;
 553           outval4 += masque[ base +3 ]*pix ;
 554           outval5 += masque[ base +2 ]*pix ;
 555           outval6 += masque[ base +1 ]*pix ;
 556           outval7 += masque[ base    ]*pix ;
 557
 558           pix = tex2D(tex_img_inc, j+1, i-2+ic) ;
 559           outval0 += masque[ base +3 ]*pix ;
 560           outval1 += masque[ base +2 ]*pix ;
 561           outval2 += masque[ base +1 ]*pix ;
 562           outval3 += masque[ base    ]*pix ;
 563
 564           pix = tex2D(tex_img_inc, j, i-2+ic) ;
 565           outval0 += masque[ base +2 ]*pix ;
 566           outval1 += masque[ base +1 ]*pix ;
 567           outval2 += masque[ base    ]*pix ;
 568
 569           pix = tex2D(tex_img_inc, j-1, i-2+ic) ;
 570           outval0 += masque[ base +1 ]*pix ;
 571           outval1 += masque[ base    ]*pix ;
 572
 573           pix = tex2D(tex_img_inc, j-2, i-2+ic) ;
 574           outval0 += masque[ base    ]*pix ;
 575
 576           pix = tex2D(tex_img_inc, j+6, i-2+ic) ;
 577           outval7 += masque[ base +1 ]*pix ;
 578           outval6 += masque[ base +2 ]*pix ;
 579           outval5 += masque[ base +3 ]*pix ;
 580           outval4 += masque[ base +4 ]*pix ;
 581
 582           pix = tex2D(tex_img_inc, j+7, i-2+ic) ;
 583           outval7 += masque[ base +2 ]*pix ;
 584           outval6 += masque[ base +3 ]*pix ;
 585           outval5 += masque[ base +4 ]*pix ;
 586
 587           pix = tex2D(tex_img_inc, j+8, i-2+ic) ;
 588           outval7 += masque[ base +3 ]*pix ;
 589           outval6 += masque[ base +4 ]*pix ;
 590
 591           pix = tex2D(tex_img_inc, j+9, i-2+ic) ;
 592           outval7 += masque[ base +4 ]*pix ;
 593         }
 594
 595   // les 8 pixels
 596   base = __umul24(i, j_dim) + j ;
 597   output[ base++ ] = outval0 ;
 598   output[ base++ ] = outval1 ;
 599   output[ base++ ] = outval2 ;
 600   output[ base++ ] = outval3;
 601   output[ base++ ] = outval4;
 602   output[ base++ ] = outval5;
 603   output[ base++ ] = outval6 ;
 604   output[ base++ ] = outval7 ;
 605
 606 }
 607
 608 __global__ void kernel_convoGene8x8pL7( unsigned char  *output, int j_dim )
 609 {
 610   int ic, base ;
 611   int L=7 ;
 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
 616   // coordonnees absolues du point de base en haut a gauche
 617   int j = ( __umul24( blockIdx.x, blockDim.x) + threadIdx.x)<< 3 ;
 618   int i = ( __umul24( blockIdx.y, blockDim.y) + threadIdx.y) ;
 619
 620   //zone commune
 621   //forme generique contrib pixel => masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic)
 622   for (ic=0 ; ic<L ; ic++)
 623         {
 624           base = __umul24(ic,L) ;
 625           pix = tex2D(tex_img_inc, j+3, i-3+ic) ;
 626           outval0 += masque[ base +6 ]*pix ;
 627           outval1 += masque[ base +5 ]*pix ;
 628           outval2 += masque[ base +4 ]*pix ;
 629           outval3 += masque[ base +3 ]*pix ;
 630           outval4 += masque[ base +2 ]*pix ;
 631           outval5 += masque[ base +1 ]*pix ;
 632           outval6 += masque[ base    ]*pix ;
 633           pix = tex2D(tex_img_inc, j+4, i-3+ic) ;
 634           outval1 += masque[ base +6 ]*pix ;
 635           outval2 += masque[ base +5 ]*pix ;
 636           outval3 += masque[ base +4 ]*pix ;
 637           outval4 += masque[ base +3 ]*pix ;
 638           outval5 += masque[ base +2 ]*pix ;
 639           outval6 += masque[ base +1 ]*pix ;
 640           outval7 += masque[ base    ]*pix ;
 641
 642           pix = tex2D(tex_img_inc, j+2, i-3+ic) ;
 643           outval0 += masque[ base +5 ]*pix ;
 644           outval1 += masque[ base +4 ]*pix ;
 645           outval2 += masque[ base +3 ]*pix ;
 646           outval3 += masque[ base +2 ]*pix ;
 647           outval4 += masque[ base +1 ]*pix ;
 648           outval5 += masque[ base    ]*pix ;
 649           pix = tex2D(tex_img_inc, j+1, i-3+ic) ;
 650           outval0 += masque[ base +4 ]*pix ;
 651           outval1 += masque[ base +3 ]*pix ;
 652           outval2 += masque[ base +2 ]*pix ;
 653           outval3 += masque[ base +1 ]*pix ;
 654           outval4 += masque[ base    ]*pix ;
 655           pix = tex2D(tex_img_inc, j, i-3+ic) ;
 656           outval0 += masque[ base +3 ]*pix ;
 657           outval1 += masque[ base +2 ]*pix ;
 658           outval2 += masque[ base +1 ]*pix ;
 659           outval3 += masque[ base    ]*pix ;
 660           pix = tex2D(tex_img_inc, j-1, i-3+ic) ;
 661           outval0 += masque[ base +2 ]*pix ;
 662           outval1 += masque[ base +1 ]*pix ;
 663           outval2 += masque[ base    ]*pix ;
 664           pix = tex2D(tex_img_inc, j-2, i-3+ic) ;
 665           outval0 += masque[ base +1 ]*pix ;
 666           outval1 += masque[ base    ]*pix ;
 667           pix = tex2D(tex_img_inc, j-3, i-3+ic) ;
 668           outval0 += masque[ base ]*pix ;
 669
 670
 671           pix = tex2D(tex_img_inc, j+5, i-3+ic) ;
 672           outval7 += masque[ base +1 ]*pix ;
 673           outval6 += masque[ base +2 ]*pix ;
 674           outval5 += masque[ base +3 ]*pix ;
 675           outval4 += masque[ base +4 ]*pix ;
 676           outval3 += masque[ base +5 ]*pix ;
 677           outval2 += masque[ base +6 ]*pix ;
 678           pix = tex2D(tex_img_inc, j+6, i-3+ic) ;
 679           outval7 += masque[ base +2 ]*pix ;
 680           outval6 += masque[ base +3 ]*pix ;
 681           outval5 += masque[ base +4 ]*pix ;
 682           outval4 += masque[ base +5 ]*pix ;
 683           outval3 += masque[ base +6 ]*pix ;
 684           pix = tex2D(tex_img_inc, j+7, i-3+ic) ;
 685           outval7 += masque[ base +3 ]*pix ;
 686           outval6 += masque[ base +4 ]*pix ;
 687           outval5 += masque[ base +5 ]*pix ;
 688           outval4 += masque[ base +6 ]*pix ;
 689           pix = tex2D(tex_img_inc, j+8, i-3+ic) ;
 690           outval7 += masque[ base +4 ]*pix ;
 691           outval6 += masque[ base +5 ]*pix ;
 692           outval5 += masque[ base +6 ]*pix ;
 693           pix = tex2D(tex_img_inc, j+9, i-3+ic) ;
 694           outval7 += masque[ base +5 ]*pix ;
 695           outval6 += masque[ base +6 ]*pix ;
 696           pix = tex2D(tex_img_inc, j+10, i-3+ic) ;
 697           outval7 += masque[ base +6 ]*pix ;
 698
 699         }
 700
 701   // les 8 pixels
 702   base = __umul24(i, j_dim) + j ;
 703   output[ base++ ] = outval0 ;
 704   output[ base++ ] = outval1 ;
 705   output[ base++ ] = outval2 ;
 706   output[ base++ ] = outval3;
 707   output[ base++ ] = outval4;
 708   output[ base++ ] = outval5;
 709   output[ base++ ] = outval6 ;
 710   output[ base++ ] = outval7 ;
 711 }
 712
 713
 714
 715 // convolution non separable 4 pixels/thread
 716 // image en texture et noyau en mem constante
 717 __global__ void kernel_convoGene8rx16p( unsigned char  *output, int j_dim, int r )
 718 {
 719   int ic, jc ;
 720   int L=2*r+1, baseIdx ;
 721   unsigned char pix ;
 722   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
 723   float outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
 724   float outval8=0.0, outval9=0.0, outval10=0.0, outval11=0.0 ;
 725   float outval12=0.0, outval13=0.0, outval14=0.0, outval15=0.0 ;
 726
 727   // coordonnees absolues du point de base en haut a gauche
 728   int j = ( __umul24( blockIdx.x, blockDim.x) + threadIdx.x)<< 3 ;
 729   int i = ( __umul24( blockIdx.y, blockDim.y) + threadIdx.y)<< 1 ;
 730
 731   //zone commune
 732   //forme generique contrib pixel => masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic)
 733   for (ic=1 ; ic<L ; ic++)
 734         {
 735           for(jc=1 ; jc<L ; jc++)
 736                 {
 737                   pix = tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
 738                   baseIdx = __umul24(ic-1,L)+jc-1 ;
 739                   outval0 += masque[ baseIdx+L+1 ]*pix ;
 740                   outval1 += masque[ baseIdx+L   ]*pix ;
 741                   outval2 += masque[ baseIdx+1   ]*pix ;
 742                   outval3 += masque[ baseIdx     ]*pix ;
 743
 744                   pix = tex2D(tex_img_inc, j-r+jc+2, i-r+ic) ;
 745                   outval4 += masque[ baseIdx+L+1  ]*pix ;
 746                   outval5 += masque[ baseIdx+L    ]*pix ;
 747                   outval6 += masque[ baseIdx+1    ]*pix ;
 748                   outval7 += masque[ baseIdx      ]*pix ;
 749
 750                   pix = tex2D(tex_img_inc, j-r+jc+4, i-r+ic) ;
 751                   outval8 += masque[ baseIdx+L+1 ]*pix ;
 752                   outval9 += masque[ baseIdx+L   ]*pix ;
 753                   outval10+= masque[ baseIdx+1   ]*pix ;
 754                   outval11+= masque[ baseIdx     ]*pix ;
 755
 756                   pix = tex2D(tex_img_inc, j-r+jc+6, i-r+ic) ;
 757                   outval12 += masque[ baseIdx+L+1 ]*pix ;
 758                   outval13 += masque[ baseIdx+L   ]*pix ;
 759                   outval14 += masque[ baseIdx+1   ]*pix ;
 760                   outval15 += masque[ baseIdx     ]*pix ;
 761
 762                 }
 763         }
 764   // les blocs peripheriques verticaux
 765   for (ic=1 ; ic<L ; ic++)
 766         {
 767           outval0 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r, i-r+ic) ; // bande A pour b0
 768           outval2 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r, i-r+ic) ; // bande A pour b2
 769
 770           outval1 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+1, i-r+ic) ; //bande C pour b1
 771           outval3 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+1, i-r+ic) ; //bande C pour b3
 772
 773           outval4 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r+2, i-r+ic) ; // bande A pour b0
 774           outval6 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r+2, i-r+ic) ; // bande A pour b2
 775
 776           outval5 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+3, i-r+ic) ; //bande C pour b1
 777           outval7 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+3, i-r+ic) ; //bande C pour b3
 778
 779           outval8 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r+4, i-r+ic) ; // bande A pour b0
 780           outval10 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r+4, i-r+ic) ; // bande A pour b2
 781
 782           outval9 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+5, i-r+ic) ; //bande C pour b1
 783           outval11 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+5, i-r+ic) ; //bande C pour b3
 784
 785           outval12 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r+6, i-r+ic) ; // bande A pour b0
 786           outval14 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r+6, i-r+ic) ; // bande A pour b2
 787
 788           outval13 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+7, i-r+ic) ; //bande C pour b1
 789           outval15 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+7, i-r+ic) ; //bande C pour b3
 790         }
 791   // les blocs peripheriques horizontaux
 792   for (jc=1 ; jc<L ; jc++)
 793         {
 794           outval0 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc, i-r) ; // bande B pour b0
 795           outval1 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc, i-r) ; // bande B pour b1
 796
 797           outval2 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc, i+r+1) ; //bande D pour b2
 798           outval3 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc, i+r+1) ; //bande D pour b3
 799
 800           outval4 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc+2, i-r) ; // bande B pour b0
 801           outval5 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc+2, i-r) ; // bande B pour b1
 802
 803           outval6 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc+2, i+r+1) ; //bande D pour b2
 804           outval7 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc+2, i+r+1) ; //bande D pour b3
 805
 806           outval8 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc+4, i-r) ; // bande B pour b0
 807           outval9 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc+4, i-r) ; // bande B pour b1
 808
 809           outval10 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc+4, i+r+1) ; //bande D pour b2
 810           outval11 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc+4, i+r+1) ; //bande D pour b3
 811
 812           outval12 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc+6, i-r) ; // bande B pour b0
 813           outval13 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc+6, i-r) ; // bande B pour b1
 814
 815           outval14 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc+6, i+r+1) ; //bande D pour b2
 816           outval15 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc+6, i+r+1) ; //bande D pour b3
 817         }
 818
 819   //les coins
 820   {
 821         outval0 += masque[ 0                   ]*tex2D(tex_img_inc, j-r  , i-r) ;   // haut gauche pour  b0
 822         outval1 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+1, i-r) ;   // haut droit pour b1
 823         outval2 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r  , i+r+1) ; // bas gauche pour b2
 824         outval3 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+1, i+r+1) ; //bas droit pour b3
 825
 826         outval4 += masque[ 0                   ]*tex2D(tex_img_inc, j-r+2, i-r) ;   // haut gauche pour  b0
 827         outval5 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+3, i-r) ;   // haut droit pour b1
 828         outval6 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r+2, i+r+1) ; // bas gauche pour b2
 829         outval7 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+3, i+r+1) ; //bas droit pour b3
 830
 831         outval8 += masque[ 0                   ]*tex2D(tex_img_inc, j-r+4, i-r) ;   // haut gauche pour  b0
 832         outval9 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+5, i-r) ;   // haut droit pour b1
 833         outval10 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r+4, i+r+1) ; // bas gauche pour b2
 834         outval11 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+5, i+r+1) ; //bas droit pour b3
 835
 836         outval12 += masque[ 0                   ]*tex2D(tex_img_inc, j-r+6, i-r) ;   // haut gauche pour  b0
 837         outval13 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+7, i-r) ;   // haut droit pour b1
 838         outval14 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r+6, i+r+1) ; // bas gauche pour b2
 839         outval15 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+7, i+r+1) ; //bas droit pour b3
 840   }
 841
 842
 843   // les 2 pixels hauts
 844   int base = __umul24(i, j_dim) +j ;
 845   output[ base++   ] = outval0 ;
 846   output[ base++ ] = outval1 ;
 847   output[ base++ ] = outval4 ;
 848   output[ base++ ] = outval5 ;
 849   output[ base++ ] = outval8 ;
 850   output[ base++ ] = outval9 ;
 851   output[ base++ ] = outval12 ;
 852   output[ base++ ] = outval13 ;
 853   // les 2 pixels bas
 854   base = __umul24(i+1, j_dim) +j ;
 855   output[ base++ ] = outval2 ;
 856   output[ base++ ] = outval3 ;
 857   output[ base++ ] = outval6 ;
 858   output[ base++ ] = outval7 ;
 859   output[ base++ ] = outval10 ;
 860   output[ base++ ] = outval11 ;
 861   output[ base++ ] = outval14 ;
 862   output[ base++ ] = outval15 ;
 863
 864 }
 865
 866
 867
 868 // convolution non separable
 869 // image en texture et noyau en mem constante
 870 __global__ void kernel_convoGene8r( unsigned char  *output, int i_dim, int j_dim, int r)
 871 {
 872   int ic, jc ;
 873   int L=2*r+1 ;
 874   float outval0=0.0 ;
 875
 876   // coordonnees absolues du point de base
 877   int j = __mul24( blockIdx.x,blockDim.x ) + threadIdx.x ;
 878   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
 879
 880   for (ic=-r ; ic<=r ; ic++)
 881         for(jc=-r ; jc<=r ; jc++)
 882           outval0 += masque[ __umul24(ic,L)+jc+r ]*tex2D(tex_img_inc, j+jc, i+ic) ;
 883
 884
 885   output[ __umul24(i, j_dim) + j ] = outval0 ;
 886 }
 887
 888 // convolution non separable
 889 // image en texture et noyau en mem constante
 890 __global__ void kernel_convoGeneTex8r( unsigned char  *output, int i_dim, int j_dim, int r)
 891 {
 892   int ic, jc ;
 893   int L=2*r+1 ;
 894   float outval0=0.0 ;
 895
 896   // coordonnees absolues du point de base
 897   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
 898   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
 899
 900   for (ic=0 ; ic<L ; ic++)
 901         for(jc=0 ; jc<L ; jc++)
 902           outval0 += tex1D(tex_noyau, __mul24(ic,L)+jc )*tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
 903
 904   output[ __mul24(i, j_dim) + j ] = outval0 ;
 905 }
 906
 907
 908 __global__ void kernel_convoGene3Reg8( unsigned char  *output, int i_dim, int j_dim )
 909 {
 910   float outval0=0.0 ;
 911   float n0,n1,n2,n3,n4,n5,n6,n7,n8 ;
 912
 913   n0 = (1.0/45) ;
 914   n1 = (2.0/45) ;
 915   n2 = (3.0/45) ;
 916   n3 = (4.0/45) ;
 917   n4 = (5.0/45) ;
 918   n5 = (6.0/45) ;
 919   n6 = (7.0/45) ;
 920   n7 = (8.0/45) ;
 921   n8 = (9.0/45) ;
 922
 923   // coordonnees absolues du point de base
 924   int j = __mul24(blockIdx.x, blockDim.x) + threadIdx.x ;
 925   int i = __mul24(blockIdx.y, blockDim.y) + threadIdx.y ;
 926
 927   outval0 = n8*tex2D(tex_img_inc, j-1, i-1 )
 928               + n7*tex2D(tex_img_inc, j  , i-1 )
 929               + n6*tex2D(tex_img_inc, j+1, i-1 )
 930               + n5*tex2D(tex_img_inc, j-1, i   )
 931               + n4*tex2D(tex_img_inc, j  , i   )
 932           + n3*tex2D(tex_img_inc, j+1, i   )
 933               + n2*tex2D(tex_img_inc, j-1, i+1 )
 934               + n1*tex2D(tex_img_inc, j  , i+1 )
 935               + n0*tex2D(tex_img_inc, j+1, i+1 ) ;
 936
 937   output[ __mul24(i, j_dim) + j  ] = (unsigned char) outval0 ;
 938 }
 939
 940
 941 // convolution non separable
 942 // image en texture et noyau en registres
 943 __global__ void kernel_convoGene5Reg8( unsigned char  *output, int i_dim, int j_dim)
 944 {
 945   float outval0=0.0 ;
 946   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 ;
 947
 948   n0 = (1.0/25) ;
 949   n1 = (1.0/25) ;
 950   n2 = (1.0/25) ;
 951   n3 = (1.0/25) ;
 952   n4 = (1.0/25) ;
 953   n5 = (1.0/25) ;
 954   n6 = (1.0/25) ;
 955   n7 = (1.0/25) ;
 956   n8 = (1.0/25) ;
 957   n9 = (1.0/25) ;
 958   n10 = (1.0/25) ;
 959   n11 = (1.0/25) ;
 960   n12 = (1.0/25) ;
 961   n13 = (1.0/25) ;
 962   n14 = (1.0/25) ;
 963   n15 = (1.0/25) ;
 964   n16 = (1.0/25) ;
 965   n17 = (1.0/25) ;
 966   n18 = (1.0/25) ;
 967   n19 = (1.0/25) ;
 968   n20 = (1.0/25) ;
 969   n21 = (1.0/25) ;
 970   n22 = (1.0/25) ;
 971   n23 = (1.0/25) ;
 972   n24 = (1.0/25) ;
 973
 974   // coordonnees absolues du point de base
 975   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
 976   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
 977
 978   outval0 = n0*tex2D(tex_img_inc, j-2, i-2 )
 979         + n1*tex2D(tex_img_inc, j-1, i-2 )
 980         + n2*tex2D(tex_img_inc, j  , i-2 )
 981         + n3*tex2D(tex_img_inc, j+1, i-2 )
 982         + n4*tex2D(tex_img_inc, j+2, i-2 )
 983     + n5*tex2D(tex_img_inc, j-2, i-1 )
 984         + n6*tex2D(tex_img_inc, j-1, i-1 )
 985         + n7*tex2D(tex_img_inc, j  , i-1 )
 986         + n8*tex2D(tex_img_inc, j+1, i-1 )
 987         + n9*tex2D(tex_img_inc, j+2, i-1 )
 988   + n10*tex2D(tex_img_inc, j-2, i )
 989   + n11*tex2D(tex_img_inc, j-1, i )
 990   + n12*tex2D(tex_img_inc, j  , i )
 991   + n13*tex2D(tex_img_inc, j+1, i )
 992   + n14*tex2D(tex_img_inc, j+2, i )
 993         + n15*tex2D(tex_img_inc, j-2, i+1 )
 994         + n16*tex2D(tex_img_inc, j-1, i+1 )
 995         + n17*tex2D(tex_img_inc, j  , i+1 )
 996         + n18*tex2D(tex_img_inc, j+1, i+1 )
 997         + n19*tex2D(tex_img_inc, j+2, i+1 )
 998         + n20*tex2D(tex_img_inc, j-2, i+2 )
 999   + n21*tex2D(tex_img_inc, j-1, i+2 )
1000   + n22*tex2D(tex_img_inc, j  , i+2 )
1001   + n23*tex2D(tex_img_inc, j+1, i+2 )
1002         + n24*tex2D(tex_img_inc, j+2, i+2 ) ;
1003
1004   output[ __mul24(i, j_dim) + j  ] = outval0 ;
1005 }
1006
1007 // convolution non separable
1008 // image en texture et noyau en registres
1009 __global__ void kernel_convoGene7Reg8( unsigned char  *output, int i_dim, int j_dim)
1010 {
1011   float outval0=0.0 ;
1012   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 ;
1013   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 ;
1014
1015   n0 = (1.0/49) ;
1016   n1 = (1.0/49) ;
1017   n2 = (1.0/49) ;
1018   n3 = (1.0/49) ;
1019   n4 = (1.0/49) ;
1020   n5 = (1.0/49) ;
1021   n6 = (1.0/49) ;
1022   n7 = (1.0/49) ;
1023   n8 = (1.0/49) ;
1024   n9 = (1.0/49) ;
1025   n10 = (1.0/49) ;
1026   n11 = (1.0/49) ;
1027   n12 = (1.0/49) ;
1028   n13 = (1.0/49) ;
1029   n14 = (1.0/49) ;
1030   n15 = (1.0/49) ;
1031   n16 = (1.0/49) ;
1032   n17 = (1.0/49) ;
1033   n18 = (1.0/49) ;
1034   n19 = (1.0/49) ;
1035   n20 = (1.0/49) ;
1036   n21 = (1.0/49) ;
1037   n22 = (1.0/49) ;
1038   n23 = (1.0/49) ;
1039   n24 = (1.0/49) ;
1040
1041   n25 = (1.0/49) ;
1042   n26 = (1.0/49) ;
1043   n27 = (1.0/49) ;
1044   n28 = (1.0/49) ;
1045   n29 = (1.0/49) ;
1046   n30 = (1.0/49) ;
1047   n31 = (1.0/49) ;
1048   n32 = (1.0/49) ;
1049   n33 = (1.0/49) ;
1050   n34 = (1.0/49) ;
1051   n35 = (1.0/49) ;
1052   n36 = (1.0/49) ;
1053   n37 = (1.0/49) ;
1054   n38 = (1.0/49) ;
1055   n39 = (1.0/49) ;
1056   n40 = (1.0/49) ;
1057   n41 = (1.0/49) ;
1058   n42 = (1.0/49) ;
1059   n43 = (1.0/49) ;
1060   n44 = (1.0/49) ;
1061   n45 = (1.0/49) ;
1062   n46 = (1.0/49) ;
1063   n47 = (1.0/49) ;
1064   n48 = (1.0/49) ;
1065
1066
1067
1068   // coordonnees absolues du point de base
1069   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
1070   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1071
1072   outval0 =
1073           n0*tex2D(tex_img_inc, j-3, i-3 )
1074         + n1*tex2D(tex_img_inc, j-2, i-3 )
1075         + n2*tex2D(tex_img_inc, j-1, i-3 )
1076         + n3*tex2D(tex_img_inc, j  , i-3 )
1077         + n4*tex2D(tex_img_inc, j+1, i-3 )
1078     + n5*tex2D(tex_img_inc, j+2, i-3 )
1079         + n6*tex2D(tex_img_inc, j+3, i-3 )
1080         +  n7*tex2D(tex_img_inc, j-3, i-2 )
1081         +  n8*tex2D(tex_img_inc, j-2, i-2 )
1082         +  n9*tex2D(tex_img_inc, j-1, i-2 )
1083         + n10*tex2D(tex_img_inc, j  , i-2 )
1084         + n11*tex2D(tex_img_inc, j+1, i-2 )
1085     + n12*tex2D(tex_img_inc, j+2, i-2 )
1086         + n13*tex2D(tex_img_inc, j+3, i-2 )
1087         + n14*tex2D(tex_img_inc, j-3, i-1 )
1088         + n15*tex2D(tex_img_inc, j-2, i-1 )
1089         + n16*tex2D(tex_img_inc, j-1, i-1 )
1090         + n17*tex2D(tex_img_inc, j  , i-1 )
1091         + n18*tex2D(tex_img_inc, j+1, i-1 )
1092     + n19*tex2D(tex_img_inc, j+2, i-1 )
1093         + n20*tex2D(tex_img_inc, j+3, i-1 )
1094         + n21*tex2D(tex_img_inc, j-3, i )
1095         + n22*tex2D(tex_img_inc, j-2, i )
1096         + n23*tex2D(tex_img_inc, j-1, i )
1097         + n24*tex2D(tex_img_inc, j  , i )
1098         + n25*tex2D(tex_img_inc, j+1, i )
1099     + n26*tex2D(tex_img_inc, j+2, i )
1100         + n27*tex2D(tex_img_inc, j+3, i )
1101         + n28*tex2D(tex_img_inc, j-3, i+1 )
1102         + n29*tex2D(tex_img_inc, j-2, i+1 )
1103         + n30*tex2D(tex_img_inc, j-1, i+1 )
1104         + n31*tex2D(tex_img_inc, j  , i+1 )
1105         + n32*tex2D(tex_img_inc, j+1, i+1 )
1106     + n33*tex2D(tex_img_inc, j+2, i+1 )
1107         + n34*tex2D(tex_img_inc, j+3, i+1 )
1108         + n35*tex2D(tex_img_inc, j-3, i+2 )
1109         + n36*tex2D(tex_img_inc, j-2, i+2 )
1110         + n37*tex2D(tex_img_inc, j-1, i+2 )
1111         + n38*tex2D(tex_img_inc, j  , i+2 )
1112         + n39*tex2D(tex_img_inc, j+1, i+2 )
1113     + n40*tex2D(tex_img_inc, j+2, i+2 )
1114         + n41*tex2D(tex_img_inc, j+3, i+2 )
1115         + n42*tex2D(tex_img_inc, j-3, i+3 )
1116         + n43*tex2D(tex_img_inc, j-2, i+3 )
1117         + n44*tex2D(tex_img_inc, j-1, i+3 )
1118         + n45*tex2D(tex_img_inc, j  , i+3 )
1119         + n46*tex2D(tex_img_inc, j+1, i+3 )
1120     + n47*tex2D(tex_img_inc, j+2, i+3 )
1121         + n48*tex2D(tex_img_inc, j+3, i+3 ) ;
1122
1123
1124   output[ __mul24(i, j_dim) + j  ] = outval0 ;
1125 }
1126
1127
1128 __global__ void kernel_convoNonSep16( unsigned short*output, int i_dim, int j_dim)
1129 {
1130   int idb, ic, jc ;
1131   int r = rnoyau ;
1132   int L = 2*r+1 ;
1133   int N = L*L ;
1134   float outval0=0 ;
1135
1136   // coordonnees absolues du point de base
1137   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
1138   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1139
1140   for (idb=0 ; idb< N ; idb++)
1141         {
1142           ic = i-r + idb/L ;
1143           jc = j-r + idb - (idb/L)*L ;
1144           outval0 += ( noyau[ idb ]*tex2D(tex_img_in, jc, ic)) ;
1145         }
1146
1147   output[ __mul24(i, j_dim) + j  ] = outval0 ;
1148 }
1149
1150 // convolution non separable
1151 // image en texture et noyau en mem constante
1152 // prefetch des pixels en smem
1153 // 1 pixel par thread
1154 __global__ void kernel_convoNonSepSh( unsigned short*output, int i_dim, int j_dim)
1155 {
1156   int idb, ic, jc ;
1157   int L = 2*rnoyau+1 ;
1158   float outval0=0 ;
1159
1160   // coordonnees absolues du point de base
1161   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
1162   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1163   int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
1164
1165   // chargement en smem
1166   int idrow = threadIdx.y*(blockDim.x+2*rnoyau) ;
1167
1168   extern __shared__ int roi[];
1169
1170   // bloc 0 (en haut \E0 gauche)
1171   roi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-rnoyau, i-rnoyau) ;
1172   // bloc 1 (en haut \E0 droite)...
1173   if ( threadIdx.x < L-1 ) //...ou plutot ce qu'il en manque
1174         roi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-rnoyau, i-rnoyau ) ;
1175   // bloc 2 ( en bas \E0 gauche)
1176   if ( threadIdx.y < L-1 )
1177         {
1178           idrow = (threadIdx.y+blockDim.y)*(blockDim.x+L-1) ;
1179           roi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-rnoyau, i+blockDim.y-rnoyau ) ;
1180           //bloc 4 ( en bas \E0 doite )...
1181           if ( threadIdx.x < L-1 ) //...ou ce qu'il en manque
1182                 roi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-rnoyau, i+blockDim.y-rnoyau ) ;
1183         }
1184   __syncthreads();
1185
1186   // calculs de convolution
1187   for (ic=0 ; ic<L ; ic++)
1188         for( jc=0 ; jc<L ; jc++)
1189           outval0 += ( noyau[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(blockDim.x+L-1)) + jc+threadIdx.x]) ;
1190
1191   // 1 pixel par thread --> global mem
1192   output[ idx  ] = outval0 ;
1193
1194 }
1195
1196 __global__ void kernel_convoNonSepSh_2p(unsigned short *output, int i_dim, int j_dim)
1197 {
1198   int idb, ic, jc ;
1199   int L = 2*rnoyau+1 ;
1200   float outval0=0.0, outval1=0.0 ;
1201   int bdimX = blockDim.x<<1 ;
1202   int tidX = threadIdx.x<<1 ;
1203
1204   // coordonnees absolues du point de base
1205   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
1206   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1207   int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
1208
1209   // chargement en smem
1210   int idrow = threadIdx.y*(bdimX+L-1) ;
1211
1212   extern __shared__ int roi[];
1213
1214   // bloc 0 (en haut \E0 gauche)
1215   roi[ idrow  + tidX ] = tex2D(tex_img_ins, j-rnoyau, i-rnoyau) ;
1216   roi[ idrow  + tidX +1] = tex2D(tex_img_ins, j-rnoyau+1, i-rnoyau) ;
1217
1218   // bloc 1 (en haut \E0 droite)...
1219   if ( threadIdx.x < rnoyau ) //...ou plutot ce qu'il en manque
1220         {
1221           roi[ idrow + tidX+bdimX   ] = tex2D( tex_img_ins, j+bdimX-rnoyau, i-rnoyau ) ;
1222           roi[ idrow + tidX+bdimX +1] = tex2D( tex_img_ins, j+bdimX-rnoyau+1, i-rnoyau ) ;
1223         }
1224
1225   // bloc 2 ( en bas \E0 gauche)
1226   if ( threadIdx.y < L-1 )
1227         {
1228           idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
1229           roi[ idrow + tidX ] = tex2D( tex_img_ins, j-rnoyau, i+blockDim.y-rnoyau ) ;
1230           roi[ idrow + tidX +1] = tex2D( tex_img_ins, j-rnoyau+1, i+blockDim.y-rnoyau ) ;
1231           //bloc 4 ( en bas \E0 doite )...
1232
1233           if ( 2*threadIdx.x < L-1 ) //...ou ce qu'il en manque
1234                 {
1235                   roi[ idrow + tidX+bdimX ] = tex2D( tex_img_ins, j+bdimX-rnoyau, i+blockDim.y-rnoyau ) ;
1236                   roi[ idrow + tidX+bdimX +1] = tex2D( tex_img_ins, j+bdimX-rnoyau+1, i+blockDim.y-rnoyau ) ;
1237                   }
1238         }
1239   __syncthreads();
1240
1241   // calculs de convolution
1242   for (ic=0 ; ic<L ; ic++)
1243         for( jc=0 ; jc<L ; jc++)
1244           {
1245                 outval0 += (noyau[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(2*blockDim.x+L-1)) + jc+2*threadIdx.x]) ;
1246                 outval1 += (noyau[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(2*blockDim.x+L-1)) + jc+2*threadIdx.x+1]) ;
1247           }
1248
1249   // 1 pixel par thread --> global mem
1250   output[ idx  ] = outval0 ;
1251   output[ idx+1 ] = outval1 ;
1252 }
1253
1254 __global__ void kernel_convoNonSepSh_2p(unsigned char *output, int i_dim, int j_dim, int r)
1255 {
1256   int idb, ic, jc ;
1257   int L = 2*r+1 ;
1258   float outval0=0.0, outval1=0.0 ;
1259   int bdimX = blockDim.x<<1 ;
1260   int tidX = threadIdx.x<<1 ;
1261
1262   // coordonnees absolues du point de base
1263   int j = __umul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
1264   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1265   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
1266
1267   // chargement en smem
1268   int idrow = threadIdx.y*(bdimX+L-1) ;
1269
1270   extern __shared__ int roi[];
1271
1272   // bloc 0 (en haut \E0 gauche)
1273   roi[ idrow  + tidX    ] = tex2D( tex_img_inc, j-r  , i-r) ;
1274   roi[ idrow  + tidX +1 ] = tex2D( tex_img_inc, j-r+1, i-r) ;
1275
1276   // bloc 1 (en haut \E0 droite)...
1277   if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
1278         {
1279           roi[ idrow + bdimX + tidX    ] = tex2D( tex_img_inc, j+bdimX-r  , i-r ) ;
1280           roi[ idrow + bdimX + tidX +1 ] = tex2D( tex_img_inc, j+bdimX-r+1, i-r ) ;
1281           }
1282
1283   // bloc 2 ( en bas \E0 gauche)
1284   if ( threadIdx.y < L-1 )
1285         {
1286           idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
1287           roi[ idrow + tidX ] = tex2D( tex_img_inc, j-r, i+blockDim.y-r ) ;
1288           roi[ idrow + tidX +1] = tex2D( tex_img_inc, j-r+1, i+blockDim.y-r ) ;
1289           //bloc 4 ( en bas \E0 doite )...
1290
1291           if ( threadIdx.x < r ) //...ou ce qu'il en manque
1292                 {
1293                   roi[ idrow + tidX+bdimX   ] = tex2D( tex_img_inc, j+bdimX-r  , i+blockDim.y-r ) ;
1294                   roi[ idrow + tidX+bdimX +1] = tex2D( tex_img_inc, j+bdimX-r+1, i+blockDim.y-r ) ;
1295                   }
1296         }
1297   __syncthreads();
1298
1299   // calculs de convolution
1300   for (ic=0 ; ic<L ; ic++)
1301         for( jc=0 ; jc<L ; jc++)
1302           {
1303                 outval0 += (masque[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX]) ;
1304                 outval1 += (masque[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX+1]) ;
1305           }
1306
1307   // 1 pixel par thread --> global mem
1308   output[ idx  ] = outval0 ;
1309   output[ idx+1 ] = outval1 ;
1310 }
1311
1312 //4 pixels en ligne par thread
1313 __global__ void kernel_convoNonSepSh_4p(unsigned char *output, int i_dim, int j_dim, int r)
1314 {
1315   int idb, ic, jc ;
1316   int L = 2*r+1 ;
1317   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
1318   int bdimX = blockDim.x<<2 ;
1319   int tidX = threadIdx.x<<2 ;
1320
1321   // coordonnees absolues du point de base
1322   int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<2 ;
1323   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1324   int j0= __umul24(blockIdx.x,blockDim.x)<<2 ;
1325
1326   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
1327
1328
1329   // chargement en smem
1330   int idrow = threadIdx.y*(bdimX+L-1) ;
1331
1332   extern __shared__ unsigned char roi4p[];
1333
1334   // bloc 0 (en haut \E0 gauche)
1335   roi4p[ idrow  + tidX   ] = tex2D(tex_img_inc, j-r  , i-r) ;
1336   roi4p[ idrow  + tidX +1] = tex2D(tex_img_inc, j-r+1, i-r) ;
1337   roi4p[ idrow  + tidX +2] = tex2D(tex_img_inc, j-r+2, i-r) ;
1338   roi4p[ idrow  + tidX +3] = tex2D(tex_img_inc, j-r+3, i-r) ;
1339
1340   // bloc 1 (en haut \E0 droite)...
1341   if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
1342         {
1343           roi4p[ idrow + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r ) ;
1344           roi4p[ idrow + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r ) ;
1345         }
1346   // bloc 2 ( en bas \E0 gauche)
1347   if ( threadIdx.y < L-1 )
1348         {
1349           idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
1350           roi4p[ idrow + tidX   ] = tex2D( tex_img_inc, j-r  , i+blockDim.y-r ) ;
1351           roi4p[ idrow + tidX +1] = tex2D( tex_img_inc, j-r+1, i+blockDim.y-r ) ;
1352           roi4p[ idrow + tidX +2] = tex2D( tex_img_inc, j-r+2, i+blockDim.y-r ) ;
1353           roi4p[ idrow + tidX +3] = tex2D( tex_img_inc, j-r+3, i+blockDim.y-r ) ;
1354
1355           //bloc 4 ( en bas \E0 doite )...
1356           if ( threadIdx.x < r ) //...ou ce qu'il en manque
1357                 {
1358                   roi4p[ idrow + bdimX +threadIdx.x   ] = tex2D( tex_img_inc, j0-r +bdimX +threadIdx.x, i+blockDim.y-r ) ;
1359                   roi4p[ idrow + bdimX +threadIdx.x +r] = tex2D( tex_img_inc, j0   +bdimX +threadIdx.x, i+blockDim.y-r ) ;
1360                 }
1361         }
1362   __syncthreads();
1363
1364   // calculs de convolution
1365   for (ic=0 ; ic<L ; ic++)
1366         for( jc=0 ; jc<L ; jc++)
1367           {
1368                 int baseRoi = __umul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX ;
1369                 float valMask = masque[ __umul24(ic,L)+jc ] ;
1370
1371                 outval0 += valMask*roi4p[ baseRoi ] ;
1372                 outval1 += valMask*roi4p[ baseRoi +1 ] ;
1373                 outval2 += valMask*roi4p[ baseRoi +2 ] ;
1374                 outval3 += valMask*roi4p[ baseRoi +3 ] ;
1375           }
1376
1377   // 1 pixel par thread --> global mem
1378   output[ idx   ] = outval0 ;
1379   output[ idx+1 ] = outval1 ;
1380   output[ idx+2 ] = outval2 ;
1381   output[ idx+3 ] = outval3 ;
1382 }
1383
1384
1385 //4 pixels en carre par thread
1386 __global__ void kernel_convoNonSepSh_4pcarre(unsigned char *output, int i_dim, int j_dim, int r)
1387 {
1388   int idb, ic, jc ;
1389   int L = 2*r+1 ;
1390   float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
1391   int bdimX = blockDim.x<<1 ;
1392   int tidX = threadIdx.x<<1 ;
1393   int bdimY = blockDim.y<<1 ;
1394   int tidY = threadIdx.y<<1 ;
1395
1396   // coordonnees absolues du point de base
1397   int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<1 ;
1398   int i = (__umul24( blockIdx.y, blockDim.y) + threadIdx.y)<<1 ;
1399   int j0= __umul24(blockIdx.x,blockDim.x)<<1 ;
1400   int i0= __umul24(blockIdx.y,blockDim.y)<<1 ;
1401
1402   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
1403
1404
1405   // chargement en smem
1406   int idrowBase   = tidY*(bdimX+L-1) ;
1407   int idrowSecond = (tidY+1)*(bdimX+L-1) ;
1408
1409
1410   extern __shared__ unsigned char roi4p2[];
1411
1412   // bloc 0 (en haut \E0 gauche)
1413   roi4p2[ idrowBase  + tidX   ] = tex2D(tex_img_inc, j-r  , i-r  ) ;
1414   roi4p2[ idrowBase  + tidX +1] = tex2D(tex_img_inc, j-r+1, i-r  ) ;
1415   roi4p2[ idrowSecond+ tidX   ] = tex2D(tex_img_inc, j-r  , i-r+1) ;
1416   roi4p2[ idrowSecond+ tidX +1] = tex2D(tex_img_inc, j-r+1, i-r+1) ;
1417
1418   // bloc 1 (en haut \E0 droite)...
1419   if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
1420         {
1421           roi4p2[ idrowBase + bdimX   + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r ) ;
1422           roi4p2[ idrowBase + bdimX   + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r ) ;
1423           roi4p2[ idrowSecond + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r+1 ) ;
1424           roi4p2[ idrowSecond + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r+1 ) ;
1425         }
1426   // bloc 2 ( en bas \E0 gauche)
1427   if ( threadIdx.y < L-1 )
1428         {
1429           idrowBase   = (bdimY + threadIdx.y)*(bdimX+L-1) ;
1430           roi4p2[ idrowBase + tidX     ] = tex2D( tex_img_inc, j-r  , i0-r +bdimY +threadIdx.y ) ;
1431           roi4p2[ idrowBase + tidX   +1] = tex2D( tex_img_inc, j-r+1, i0-r +bdimY +threadIdx.y ) ;
1432
1433           //bloc 4 ( en bas \E0 doite )...
1434           if ( threadIdx.x < r ) //...ou ce qu'il en manque
1435                 {
1436                   roi4p2[ idrowBase + bdimX + threadIdx.x    ]   = tex2D( tex_img_inc, j0-r +bdimX +threadIdx.x  , i0-r +bdimY +threadIdx.y ) ;
1437                   roi4p2[ idrowBase + bdimX + threadIdx.x +r ]   = tex2D( tex_img_inc, j0   +bdimX +threadIdx.x  , i0-r +bdimY +threadIdx.y ) ;
1438                 }
1439         }
1440   __syncthreads();
1441
1442   // calculs de convolution
1443   for (ic=0 ; ic<L ; ic++)
1444         for( jc=0 ; jc<L ; jc++)
1445           {
1446                 int baseRoi = __umul24(ic+tidY,(bdimX+L-1)) + jc+tidX ;
1447                 float valMask = masque[ __umul24(ic,L)+jc ] ;
1448
1449                 outval0 += valMask*roi4p2[ baseRoi ] ;
1450                 outval1 += valMask*roi4p2[ baseRoi +1 ] ;
1451                 outval2 += valMask*roi4p2[ baseRoi +bdimX+L-1 ] ;
1452                 outval3 += valMask*roi4p2[ baseRoi +bdimX+L ] ;
1453           }
1454
1455   // 1 pixel par thread --> global mem
1456   output[ idx   ] = outval0 ;
1457   output[ idx+1 ] = outval1 ;
1458   output[ idx+j_dim ] = outval2 ;
1459   output[ idx+j_dim+1 ] = outval3 ;
1460 }
1461
1462
1463 //8 pixels en ligne par thread
1464 __global__ void kernel_convoNonSepSh_8p(unsigned char *output, int i_dim, int j_dim, int r)
1465 {
1466   int idb, ic, jc;
1467   int L = 2*r+1 ;
1468   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 ;
1469   int bdimX = blockDim.x<<3 ;
1470   int tidX = threadIdx.x<<3 ;
1471
1472   // coordonnees absolues du point de base
1473   int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<3 ;
1474   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1475   int j0= __umul24(blockIdx.x,blockDim.x)<<3 ;
1476
1477   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
1478
1479
1480   // chargement en smem
1481   int idrow = threadIdx.y*(bdimX+L-1) ;
1482
1483   extern __shared__ unsigned char roi8p[];
1484
1485   // bloc 0 (en haut \E0 gauche)
1486   for (int p=0; p<8; p++)
1487   roi8p[ idrow  + tidX +p ] = tex2D(tex_img_inc, j-r+p  , i-r) ;
1488
1489   // bloc 1 (en haut \E0 droite)...
1490   if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
1491         {
1492           roi8p[ idrow + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r ) ;
1493           roi8p[ idrow + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r ) ;
1494         }
1495   // bloc 2 ( en bas \E0 gauche)
1496   if ( threadIdx.y < L-1 )
1497         {
1498           idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
1499           for (int p=0; p<8; p++)
1500                 roi8p[ idrow + tidX +p  ] = tex2D( tex_img_inc, j-r+p  , i+blockDim.y-r ) ;
1501
1502           //bloc 4 ( en bas \E0 doite )...
1503           if ( threadIdx.x < r ) //...ou ce qu'il en manque
1504                 {
1505                   roi8p[ idrow + bdimX +threadIdx.x   ] = tex2D( tex_img_inc, j0-r +bdimX +threadIdx.x, i+blockDim.y-r ) ;
1506                   roi8p[ idrow + bdimX +threadIdx.x +r] = tex2D( tex_img_inc, j0   +bdimX +threadIdx.x, i+blockDim.y-r ) ;
1507                 }
1508         }
1509   __syncthreads();
1510
1511   // calculs de convolution
1512   for (ic=0 ; ic<L ; ic++)
1513         for( jc=0 ; jc<L ; jc++)
1514           {
1515                 int baseRoi = __umul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX ;
1516                 float valMask = masque[ __umul24(ic,L)+jc ] ;
1517                 outval0 += valMask*roi8p[ baseRoi ] ;
1518                 outval1 += valMask*roi8p[ baseRoi +1 ] ;
1519                 outval2 += valMask*roi8p[ baseRoi +2 ] ;
1520                 outval3 += valMask*roi8p[ baseRoi +3 ] ;
1521                 outval4 += valMask*roi8p[ baseRoi +4 ] ;
1522                 outval5 += valMask*roi8p[ baseRoi +5 ] ;
1523                 outval6 += valMask*roi8p[ baseRoi +6 ] ;
1524                 outval7 += valMask*roi8p[ baseRoi +7 ] ;
1525           }
1526
1527   // 1 pixel par thread --> global mem
1528   output[ idx   ] = outval0 ;
1529   output[ idx+1 ] = outval1 ;
1530   output[ idx+2 ] = outval2 ;
1531   output[ idx+3 ] = outval3 ;
1532   output[ idx+4 ] = outval4 ;
1533   output[ idx+5 ] = outval5 ;
1534   output[ idx+6 ] = outval6 ;
1535   output[ idx+7 ] = outval7 ;
1536 }
1537
1538
1539 //16 pixels en ligne par thread
1540 __global__ void kernel_convoNonSepSh_16p(unsigned char *output, int i_dim, int j_dim, int r)
1541 {
1542   int idb, ic, jc;
1543   int L = 2*r+1 ;
1544   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 ;
1545   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 ;
1546   int bdimX = blockDim.x<<4 ;
1547   int tidX = threadIdx.x<<4 ;
1548
1549   // coordonnees absolues du point de base
1550   int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<4 ;
1551   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1552   int j0= __umul24(blockIdx.x,blockDim.x)<<4 ;
1553
1554   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
1555
1556
1557   // chargement en smem
1558   int idrow = threadIdx.y*(bdimX+L-1) ;
1559
1560   extern __shared__ unsigned char roi16p[];
1561
1562   // bloc 0 (en haut \E0 gauche)
1563   for (int p=0; p<16; p++)
1564         roi16p[ idrow  + tidX +p ] = tex2D(tex_img_inc, j-r+p  , i-r) ;
1565
1566   // bloc 1 (en haut \E0 droite)...
1567   if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
1568         {
1569           roi16p[ idrow + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r ) ;
1570           roi16p[ idrow + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r ) ;
1571         }
1572   // bloc 2 ( en bas \E0 gauche)
1573   if ( threadIdx.y < L-1 )
1574         {
1575           idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
1576           for (int p=0; p<16; p++)
1577                 roi16p[ idrow + tidX +p  ] = tex2D( tex_img_inc, j-r+p  , i+blockDim.y-r ) ;
1578
1579           //bloc 4 ( en bas \E0 doite )...
1580           if ( threadIdx.x < r ) //...ou ce qu'il en manque
1581                 {
1582                   roi16p[ idrow + bdimX +threadIdx.x   ] = tex2D( tex_img_inc, j0-r +bdimX +threadIdx.x, i+blockDim.y-r ) ;
1583                   roi16p[ idrow + bdimX +threadIdx.x +r] = tex2D( tex_img_inc, j0   +bdimX +threadIdx.x, i+blockDim.y-r ) ;
1584                 }
1585         }
1586   __syncthreads();
1587
1588   // calculs de convolution
1589   for (ic=0 ; ic<L ; ic++)
1590         for( jc=0 ; jc<L ; jc++)
1591           {
1592                 int baseRoi = __umul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX ;
1593                 float valMask = masque[ __umul24(ic,L)+jc ] ;
1594                 outval0 += valMask*roi16p[ baseRoi ] ;
1595                 outval1 += valMask*roi16p[ baseRoi +1 ] ;
1596                 outval2 += valMask*roi16p[ baseRoi +2 ] ;
1597                 outval3 += valMask*roi16p[ baseRoi +3 ] ;
1598                 outval4 += valMask*roi16p[ baseRoi +4 ] ;
1599                 outval5 += valMask*roi16p[ baseRoi +5 ] ;
1600                 outval6 += valMask*roi16p[ baseRoi +6 ] ;
1601                 outval7 += valMask*roi16p[ baseRoi +7 ] ;
1602                 outval8 += valMask*roi16p[ baseRoi +8] ;
1603                 outval9 += valMask*roi16p[ baseRoi +9 ] ;
1604                 outval10 += valMask*roi16p[ baseRoi +10 ] ;
1605                 outval11 += valMask*roi16p[ baseRoi +11 ] ;
1606                 outval12 += valMask*roi16p[ baseRoi +12 ] ;
1607                 outval13 += valMask*roi16p[ baseRoi +13 ] ;
1608                 outval14 += valMask*roi16p[ baseRoi +14 ] ;
1609                 outval15 += valMask*roi16p[ baseRoi +15 ] ;
1610           }
1611
1612   // 1 pixel par thread --> global mem
1613   output[ idx   ] = outval0 ;
1614   output[ idx+1 ] = outval1 ;
1615   output[ idx+2 ] = outval2 ;
1616   output[ idx+3 ] = outval3 ;
1617   output[ idx+4 ] = outval4 ;
1618   output[ idx+5 ] = outval5 ;
1619   output[ idx+6 ] = outval6 ;
1620   output[ idx+7 ] = outval7 ;
1621   output[ idx+8 ] = outval8 ;
1622   output[ idx+9 ] = outval9 ;
1623   output[ idx+10 ] = outval10 ;
1624   output[ idx+11 ] = outval11 ;
1625   output[ idx+12 ] = outval12 ;
1626   output[ idx+13 ] = outval13 ;
1627   output[ idx+14 ] = outval14 ;
1628   output[ idx+15 ] = outval15 ;
1629 }
1630
1631
1632
1633 // convolution non separable
1634 // image en texture et noyau en mem constante
1635 // fetch direct des pixels
1636 // 2 pixels traités par thread => meilleur débit
1637 __global__ void kernel_convoNonSep_2p( int *output, int i_dim, int j_dim)
1638 {
1639   int idb, ic, jc ;
1640   int r = rnoyau ;
1641   int L = 2*r+1 ;
1642   int N = L*L ;
1643   float outval0=0, outval1=0 ;
1644
1645   // coordonnees absolues du point de base
1646   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x, 2) ;
1647   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1648   int idx = __mul24(i, j_dim) + j ;              // indice  dans l'image
1649
1650   #pragma unroll
1651   for (idb=0 ; idb< N ; idb++)
1652         {
1653           ic = i-r + idb/L ;
1654           jc = j-r + idb - (idb/L)*L ;
1655           outval0 += ( noyau[ idb ]*tex2D(tex_img_in, jc, ic)) ;
1656           outval1 += ( noyau[ idb ]*tex2D(tex_img_in, jc+1, ic)) ;
1657         }
1658
1659   output[ idx   ] = outval0  ;
1660   output[ idx+1 ] = outval1  ;
1661
1662
1663 }
1664
1665
1666 // convolution non separable
1667 // image en texture et noyau en mem constante
1668 // fetch direct des pixels
1669 // 2 pixels traités par thread => meilleur débit
1670 __global__ void kernel_convoNonSep_2p_s( unsigned short*output, int i_dim, int j_dim)
1671 {
1672   int idb, ic, jc ;
1673   int r = rnoyau ;
1674   int L = 2*r+1 ;
1675   int N = L*L ;
1676   float outval0=0, outval1=0 ;
1677
1678   // coordonnees absolues du point de base
1679   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x, 2) ;
1680   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1681   int idx = __mul24(i, j_dim) + j ;              // indice  dans l'image
1682
1683   #pragma unroll
1684   for (idb=0 ; idb< N ; idb++)
1685         {
1686           ic = i-r + idb/L ;
1687           jc = j-r + idb - (idb/L)*L ;
1688           outval0 += ( noyau[ idb ]*tex2D(tex_img_ins, jc, ic)) ;
1689           outval1 += ( noyau[ idb ]*tex2D(tex_img_ins, jc+1, ic)) ;
1690         }
1691
1692   output[ idx   ] = outval0  ;
1693   output[ idx+1 ] = outval1  ;
1694
1695 }
1696
1697
1698
1699 // convolution separable
1700 // image en texture et noyau en mem constante
1701 __global__ void kernel_convoSep8V( unsigned char  *output, int i_dim, int j_dim, int r)
1702 {
1703   int ic ;
1704   int L=2*r+1 ;
1705   float outval0=0.0 ;
1706
1707   // coordonnees absolues du point de base
1708   int j = __mul24( blockIdx.x, blockDim.x ) + threadIdx.x ;
1709   int i = __mul24( blockIdx.y, blockDim.y ) + threadIdx.y ;
1710
1711   //vertical
1712   for (ic=0 ; ic<L ; ic++)
1713           outval0 += masque[ ic ]*tex2D(tex_img_inc, j, i-r+ic) ;
1714
1715   output[ __mul24(i, j_dim) + j ] = outval0 ;
1716 }
1717
1718
1719
1720
1721 __global__ void kernel_convoSep8H( unsigned char  *output, int i_dim, int j_dim, int r)
1722 {
1723   int jc ;
1724   int L=2*r+1 ;
1725   float outval0=0.0 ;
1726
1727   // coordonnees absolues du point de base
1728   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
1729   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
1730
1731   //horizontal
1732   for (jc=0 ; jc<L ; jc++)
1733           outval0 += masque[ L+jc ]*tex2D(tex_img_inc, j-r+jc, i) ;
1734
1735   output[ __mul24(i, j_dim) + j ] = outval0 ;
1736 }
1737
1738 // convolution separable
1739 // image en texture et noyau en mem constante
1740 __global__ void kernel_convoSep8Vx8p( unsigned char  *output, int i_dim, int j_dim, int r)
1741 {
1742   int ic, baseIdx ;
1743   int L=2*r+1 ;
1744   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 ;
1745
1746   // coordonnees absolues du point de base
1747   int j = (__mul24( blockIdx.x, blockDim.x ) + threadIdx.x)<<3 ;
1748   int i = __mul24( blockIdx.y, blockDim.y ) + threadIdx.y ;
1749
1750   //vertical
1751   for (ic=0 ; ic<L ; ic++)
1752         {
1753           float valMask = masque[ ic ] ;
1754           outval0 += valMask*tex2D(tex_img_inc, j  , i-r+ic) ;
1755           outval1 += valMask*tex2D(tex_img_inc, j+1, i-r+ic) ;
1756           outval2 += valMask*tex2D(tex_img_inc, j+2, i-r+ic) ;
1757           outval3 += valMask*tex2D(tex_img_inc, j+3, i-r+ic) ;
1758           outval4 += valMask*tex2D(tex_img_inc, j+4, i-r+ic) ;
1759           outval5 += valMask*tex2D(tex_img_inc, j+5, i-r+ic) ;
1760           outval6 += valMask*tex2D(tex_img_inc, j+6, i-r+ic) ;
1761           outval7 += valMask*tex2D(tex_img_inc, j+7, i-r+ic) ;
1762         }
1763
1764   baseIdx = __mul24(i  , j_dim) + j ;
1765   output[ baseIdx++ ] = outval0 ;
1766   output[ baseIdx++ ] = outval1 ;
1767   output[ baseIdx++ ] = outval2 ;
1768   output[ baseIdx++ ] = outval3 ;
1769   output[ baseIdx++ ] = outval4 ;
1770   output[ baseIdx++ ] = outval5 ;
1771   output[ baseIdx++ ] = outval6 ;
1772   output[ baseIdx++ ] = outval7 ;
1773 }
1774
1775
1776
1777 // convolution separable
1778 // image en texture et noyau en mem constante
1779 __global__ void kernel_convoSep8HL3x8p( unsigned char  *output, int i_dim, int j_dim)
1780 {
1781   int jc, baseIdx ;
1782   int L=3 ;
1783   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 ;
1784   float val ;
1785   // coordonnees absolues du point de base
1786   int j = (__mul24( blockIdx.x, blockDim.x ) + threadIdx.x)<<3 ;
1787   int i = __mul24( blockIdx.y, blockDim.y ) + threadIdx.y ;
1788
1789   //horizontal
1790         val = tex2D(tex_img_inc, j-1  , i) ;
1791 //p1
1792         outval0 += masque[0]*val ;
1793 //p2
1794 val = tex2D(tex_img_inc, j  , i) ;
1795         outval0 += masque[1]*val ;
1796         outval1 += masque[0]*val ;
1797 //p3
1798 val = tex2D(tex_img_inc, j+1  , i) ;
1799         outval0 += masque[2]*val ;
1800         outval1 += masque[1]*val ;
1801         outval2 += masque[0]*val ;
1802 //p4
1803 val = tex2D(tex_img_inc, j+2  , i) ;
1804         outval1 += masque[2]*val ;
1805         outval2 += masque[1]*val ;
1806         outval3 += masque[0]*val ;
1807 //p5
1808 val = tex2D(tex_img_inc, j+3  , i) ;
1809         outval2 += masque[2]*val ;
1810         outval3 += masque[1]*val ;
1811         outval4 += masque[0]*val ;
1812 //p6
1813 val = tex2D(tex_img_inc, j+4  , i) ;
1814         outval3 += masque[2]*val ;
1815         outval4 += masque[1]*val ;
1816         outval5 += masque[0]*val ;
1817 //p7
1818 val = tex2D(tex_img_inc, j+5  , i) ;
1819         outval4 += masque[2]*val ;
1820         outval5 += masque[1]*val ;
1821         outval6 += masque[0]*val ;
1822 //p8
1823 val = tex2D(tex_img_inc, j+6  , i) ;
1824         outval5 += masque[2]*val ;
1825         outval6 += masque[1]*val ;
1826         outval7 += masque[0]*val ;
1827 //p9
1828 val = tex2D(tex_img_inc, j+7  , i) ;
1829         outval6 += masque[2]*val ;
1830         outval7 += masque[1]*val ;
1831 //p10
1832 val = tex2D(tex_img_inc, j+8  , i) ;
1833         outval7 += masque[2]*val ;
1834
1835   baseIdx = __mul24(i  , j_dim) + j ;
1836   output[ baseIdx++ ] = outval0 ;
1837   output[ baseIdx++ ] = outval1 ;
1838   output[ baseIdx++ ] = outval2 ;
1839   output[ baseIdx++ ] = outval3 ;
1840   output[ baseIdx++ ] = outval4 ;
1841   output[ baseIdx++ ] = outval5 ;
1842   output[ baseIdx++ ] = outval6 ;
1843   output[ baseIdx++ ] = outval7 ;
1844 }
1845
1846
1847 // convolution separable
1848 // image en texture et noyau en mem constante
1849 __global__ void kernel_convoSep8HL3x8pG( unsigned char *input, unsigned char  *output, int i_dim, int j_dim)
1850 {
1851   int jc, baseIdx, id ;
1852   int L=3 ;
1853   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 ;
1854   float val ;
1855   // coordonnees absolues du point de base
1856   int j = (__mul24( blockIdx.x, blockDim.x ) + threadIdx.x)<<3 ;
1857   int i = __mul24( blockIdx.y, blockDim.y ) + threadIdx.y ;
1858
1859   //horizontal
1860 baseIdx = __mul24(i  , j_dim) + j ;
1861 if (baseIdx >0 ) id = baseIdx-1 ; else id = baseIdx ;
1862         val = input[id] ;
1863 //p1
1864         outval0 += masque[0]*val ;
1865 //p2
1866 val = input[baseIdx++] ;
1867         outval0 += masque[1]*val ;
1868         outval1 += masque[0]*val ;
1869 //p3
1870 val = input[baseIdx++] ;
1871         outval0 += masque[2]*val ;
1872         outval1 += masque[1]*val ;
1873         outval2 += masque[0]*val ;
1874 //p4
1875 val = input[baseIdx++] ;
1876         outval1 += masque[2]*val ;
1877         outval2 += masque[1]*val ;
1878         outval3 += masque[0]*val ;
1879 //p5
1880 val = input[baseIdx++] ;
1881         outval2 += masque[2]*val ;
1882         outval3 += masque[1]*val ;
1883         outval4 += masque[0]*val ;
1884 //p6
1885 val = input[baseIdx++] ;
1886         outval3 += masque[2]*val ;
1887         outval4 += masque[1]*val ;
1888         outval5 += masque[0]*val ;
1889 //p7
1890 val = input[baseIdx++] ;
1891         outval4 += masque[2]*val ;
1892         outval5 += masque[1]*val ;
1893         outval6 += masque[0]*val ;
1894 //p8
1895 val = input[baseIdx++] ;
1896         outval5 += masque[2]*val ;
1897         outval6 += masque[1]*val ;
1898         outval7 += masque[0]*val ;
1899 //p9
1900 val = input[baseIdx++] ;
1901         outval6 += masque[2]*val ;
1902         outval7 += masque[1]*val ;
1903 //p10
1904 val = input[baseIdx++] ;
1905         outval7 += masque[2]*val ;
1906
1907   baseIdx = __mul24(i  , j_dim) + j ;
1908   output[ baseIdx++ ] = outval0 ;
1909   output[ baseIdx++ ] = outval1 ;
1910   output[ baseIdx++ ] = outval2 ;
1911   output[ baseIdx++ ] = outval3 ;
1912   output[ baseIdx++ ] = outval4 ;
1913   output[ baseIdx++ ] = outval5 ;
1914   output[ baseIdx++ ] = outval6 ;
1915   output[ baseIdx++ ] = outval7 ;
1916 }
1917
1918 // convolution separable
1919 // image en texture et noyau en mem constante
1920 __global__ void kernel_convoSep8HL5x8pG( unsigned char *input, unsigned char  *output, int i_dim, int j_dim)
1921 {
1922   int jc, baseIdx, id ;
1923   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 ;
1924   float val ;
1925   // coordonnees absolues du point de base
1926   int j = (__mul24( blockIdx.x, blockDim.x ) + threadIdx.x)<<3 ;
1927   int i = __mul24( blockIdx.y, blockDim.y ) + threadIdx.y ;
1928
1929   //horizontal
1930 baseIdx = __mul24(i  , j_dim) + j ;
1931 if (baseIdx >1 ) id = baseIdx-2 ; else id = baseIdx ;
1932 val = input[id] ;
1933 //p1
1934         outval0 += masque[0]*val ;
1935 //p2
1936 if (baseIdx >0 ) id = baseIdx-1 ; else id = baseIdx ;
1937 val = input[id] ;
1938         outval0 += masque[1]*val ;
1939         outval1 += masque[0]*val ;
1940 //p3
1941 val = input[baseIdx++] ;
1942         outval0 += masque[2]*val ;
1943         outval1 += masque[1]*val ;
1944         outval2 += masque[0]*val ;
1945 //p4
1946 val = input[baseIdx++] ;
1947         outval0 += masque[3]*val ;
1948         outval1 += masque[2]*val ;
1949         outval2 += masque[1]*val ;
1950         outval3 += masque[0]*val ;
1951 //p5
1952 val = input[baseIdx++] ;
1953         outval0 += masque[4]*val ;
1954         outval1 += masque[3]*val ;
1955         outval2 += masque[2]*val ;
1956         outval3 += masque[1]*val ;
1957         outval4 += masque[0]*val ;
1958 //p6
1959 val = input[baseIdx++] ;
1960         outval1 += masque[4]*val ;
1961         outval2 += masque[3]*val ;
1962         outval3 += masque[2]*val ;
1963         outval4 += masque[1]*val ;
1964         outval5 += masque[0]*val ;
1965 //p7
1966 val = input[baseIdx++] ;
1967         outval2 += masque[4]*val ;
1968         outval3 += masque[3]*val ;
1969         outval4 += masque[2]*val ;
1970         outval5 += masque[1]*val ;
1971         outval6 += masque[0]*val ;
1972 //p8
1973 val = input[baseIdx++] ;
1974         outval3 += masque[4]*val ;
1975         outval4 += masque[3]*val ;
1976         outval5 += masque[2]*val ;
1977         outval6 += masque[1]*val ;
1978         outval7 += masque[0]*val ;
1979 //p9
1980 val = input[baseIdx++] ;
1981         outval4 += masque[4]*val ;
1982         outval5 += masque[3]*val ;
1983         outval6 += masque[2]*val ;
1984         outval7 += masque[1]*val ;
1985 //p10
1986 val = input[baseIdx++] ;
1987         outval5 += masque[4]*val ;
1988         outval6 += masque[3]*val ;
1989         outval7 += masque[2]*val ;
1990 //p11
1991 val = input[baseIdx++] ;
1992         outval6 += masque[4]*val ;
1993         outval7 += masque[3]*val ;
1994 //p11
1995 val = input[baseIdx++] ;
1996         outval7 += masque[4]*val ;
1997
1998
1999   baseIdx = __mul24(i  , j_dim) + j ;
2000   output[ baseIdx++ ] = outval0 ;
2001   output[ baseIdx++ ] = outval1 ;
2002   output[ baseIdx++ ] = outval2 ;
2003   output[ baseIdx++ ] = outval3 ;
2004   output[ baseIdx++ ] = outval4 ;
2005   output[ baseIdx++ ] = outval5 ;
2006   output[ baseIdx++ ] = outval6 ;
2007   output[ baseIdx++ ] = outval7 ;
2008 }
2009
2010 // convolution separable
2011 // image en texture et noyau en mem constante
2012 __global__ void kernel_convoSep8HL7x8pG( unsigned char *input, unsigned char  *output, int i_dim, int j_dim)
2013 {
2014   int jc, baseIdx, id ;
2015   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 ;
2016   float val ;
2017   // coordonnees absolues du point de base
2018   int j = (__mul24( blockIdx.x, blockDim.x ) + threadIdx.x)<<3 ;
2019   int i = __mul24( blockIdx.y, blockDim.y ) + threadIdx.y ;
2020
2021   //horizontal
2022 baseIdx = __mul24(i  , j_dim) + j ;
2023 if (baseIdx >2 ) id = baseIdx-3 ; else id = baseIdx ;
2024 val = input[id] ;
2025 //p1
2026         outval0 += masque[0]*val ;
2027 //p2
2028 if (baseIdx >1 ) id = baseIdx-2 ; else id = baseIdx ;
2029 val = input[id] ;
2030         outval0 += masque[1]*val ;
2031         outval1 += masque[0]*val ;
2032 //p3
2033 if (baseIdx >0 ) id = baseIdx-1 ; else id = baseIdx ;
2034 val = input[id] ;
2035         outval0 += masque[2]*val ;
2036         outval1 += masque[1]*val ;
2037         outval2 += masque[0]*val ;
2038 //p4
2039 val = input[baseIdx++] ;
2040         outval0 += masque[3]*val ;
2041         outval1 += masque[2]*val ;
2042         outval2 += masque[1]*val ;
2043         outval3 += masque[0]*val ;
2044 //p5
2045 val = input[baseIdx++] ;
2046         outval0 += masque[4]*val ;
2047         outval1 += masque[3]*val ;
2048         outval2 += masque[2]*val ;
2049         outval3 += masque[1]*val ;
2050         outval4 += masque[0]*val ;
2051 //p6
2052 val = input[baseIdx++] ;
2053         outval0 += masque[5]*val ;
2054         outval1 += masque[4]*val ;
2055         outval2 += masque[3]*val ;
2056         outval3 += masque[2]*val ;
2057         outval4 += masque[1]*val ;
2058         outval5 += masque[0]*val ;
2059 //p7
2060 val = input[baseIdx++] ;
2061         outval0 += masque[6]*val ;
2062         outval1 += masque[5]*val ;
2063         outval2 += masque[4]*val ;
2064         outval3 += masque[3]*val ;
2065         outval4 += masque[2]*val ;
2066         outval5 += masque[1]*val ;
2067         outval6 += masque[0]*val ;
2068 //p8
2069 val = input[baseIdx++] ;
2070         outval1 += masque[6]*val ;
2071         outval2 += masque[5]*val ;
2072         outval3 += masque[4]*val ;
2073         outval4 += masque[3]*val ;
2074         outval5 += masque[2]*val ;
2075         outval6 += masque[1]*val ;
2076         outval7 += masque[0]*val ;
2077 //p9
2078 val = input[baseIdx++] ;
2079         outval2 += masque[6]*val ;
2080         outval3 += masque[5]*val ;
2081         outval4 += masque[4]*val ;
2082         outval5 += masque[3]*val ;
2083         outval6 += masque[2]*val ;
2084         outval7 += masque[1]*val ;
2085 //p10
2086 val = input[baseIdx++] ;
2087         outval3 += masque[6]*val ;
2088         outval4 += masque[5]*val ;
2089         outval5 += masque[4]*val ;
2090         outval6 += masque[3]*val ;
2091         outval7 += masque[2]*val ;
2092 //p11
2093 val = input[baseIdx++] ;
2094         outval4 += masque[6]*val ;
2095         outval5 += masque[5]*val ;
2096         outval6 += masque[4]*val ;
2097         outval7 += masque[3]*val ;
2098 //p12
2099 val = input[baseIdx++] ;
2100         outval5 += masque[6]*val ;
2101         outval6 += masque[5]*val ;
2102         outval7 += masque[4]*val ;
2103 //p13
2104 val = input[baseIdx++] ;
2105         outval6 += masque[6]*val ;
2106         outval7 += masque[5]*val ;
2107 //p14
2108 val = input[baseIdx++] ;
2109         outval7 += masque[6]*val ;
2110
2111   baseIdx = __mul24(i  , j_dim) + j ;
2112   output[ baseIdx++ ] = outval0 ;
2113   output[ baseIdx++ ] = outval1 ;
2114   output[ baseIdx++ ] = outval2 ;
2115   output[ baseIdx++ ] = outval3 ;
2116   output[ baseIdx++ ] = outval4 ;
2117   output[ baseIdx++ ] = outval5 ;
2118   output[ baseIdx++ ] = outval6 ;
2119   output[ baseIdx++ ] = outval7 ;
2120 }
2121
2122
2123 //8 pixels en ligne par thread
2124 __global__ void kernel_convoSepShx8pV(unsigned char *output, int i_dim, int j_dim, int r)
2125 {
2126   int idb, ic, jc, p;
2127   int L = 2*r+1 ;
2128   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 ;
2129   int bdimX = blockDim.x<<3 ;
2130   int tidX = threadIdx.x<<3 ;
2131
2132   // coordonnees absolues du point de base
2133   int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<3 ;
2134   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
2135
2136   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
2137
2138
2139   // chargement en smem
2140   int idrow = threadIdx.y*bdimX ;
2141
2142   extern __shared__ unsigned char roi8p[];
2143
2144   // bloc 0 (en haut)
2145   for (p=0; p<8; p++)
2146         roi8p[ idrow  + tidX +p ] = tex2D(tex_img_inc, j+p  , i-r) ;
2147
2148
2149   // bloc 2 ( en bas)
2150   if ( threadIdx.y < L-1 )
2151         {
2152           idrow = (threadIdx.y+blockDim.y)*bdimX ;
2153           for (int p=0; p<8; p++)
2154                 roi8p[ idrow + tidX +p  ] = tex2D( tex_img_inc, j+p  , i+blockDim.y-r ) ;
2155
2156         }
2157   __syncthreads();
2158
2159   // calculs de convolution
2160   // passe verticale
2161   for (ic=0 ; ic<L ; ic++)
2162           {
2163                 int baseRoi = __umul24(ic+threadIdx.y,bdimX) + tidX ;
2164                 float valMask = masque[ ic ] ;
2165                 outval0 += valMask*roi8p[ baseRoi    ] ;
2166                 outval1 += valMask*roi8p[ baseRoi +1 ] ;
2167                 outval2 += valMask*roi8p[ baseRoi +2 ] ;
2168                 outval3 += valMask*roi8p[ baseRoi +3 ] ;
2169                 outval4 += valMask*roi8p[ baseRoi +4 ] ;
2170                 outval5 += valMask*roi8p[ baseRoi +5 ] ;
2171                 outval6 += valMask*roi8p[ baseRoi +6 ] ;
2172                 outval7 += valMask*roi8p[ baseRoi +7 ] ;
2173           }
2174
2175   // 8 pixel par thread --> global mem
2176   output[ idx   ] = outval0 ;
2177   output[ idx+1 ] = outval1 ;
2178   output[ idx+2 ] = outval2 ;
2179   output[ idx+3 ] = outval3 ;
2180   output[ idx+4 ] = outval4 ;
2181   output[ idx+5 ] = outval5 ;
2182   output[ idx+6 ] = outval6 ;
2183   output[ idx+7 ] = outval7 ;
2184 }
2185
2186 //8 pixels en ligne par thread
2187 __global__ void kernel_convoSepShx8pVG(unsigned char *input, unsigned char *output, int i_dim, int j_dim, int r)
2188 {
2189   int idb, ic, jc, p;
2190   int L = 2*r+1 ;
2191   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 ;
2192   int bdimX = blockDim.x<<3 ;
2193   int tidX = threadIdx.x<<3 ;
2194
2195   // coordonnees absolues du point de base
2196   int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<3 ;
2197   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
2198
2199   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
2200
2201
2202   // chargement en smem
2203   int idrow = threadIdx.y*bdimX ;
2204
2205   extern __shared__ unsigned char roi8p[];
2206
2207   // bloc 0 (en haut)
2208   for (p=0; p<8; p++)
2209         roi8p[ idrow  + tidX +p ] = input[ (i-r)*j_dim + j+p ] ;
2210
2211
2212   // bloc 2 ( en bas)
2213   if ( threadIdx.y < L-1 )
2214         {
2215           idrow = (threadIdx.y+blockDim.y)*bdimX ;
2216           for (int p=0; p<8; p++)
2217                 roi8p[ idrow + tidX +p  ] = input[ (i+blockDim.y-r)*j_dim + j+p ] ;
2218
2219         }
2220   __syncthreads();
2221
2222   // calculs de convolution
2223   // passe verticale
2224   for (ic=0 ; ic<L ; ic++)
2225           {
2226                 int baseRoi = __umul24(ic+threadIdx.y,bdimX) + tidX ;
2227                 float valMask = masque[ ic ] ;
2228                 outval0 += valMask*roi8p[ baseRoi    ] ;
2229                 outval1 += valMask*roi8p[ baseRoi +1 ] ;
2230                 outval2 += valMask*roi8p[ baseRoi +2 ] ;
2231                 outval3 += valMask*roi8p[ baseRoi +3 ] ;
2232                 outval4 += valMask*roi8p[ baseRoi +4 ] ;
2233                 outval5 += valMask*roi8p[ baseRoi +5 ] ;
2234                 outval6 += valMask*roi8p[ baseRoi +6 ] ;
2235                 outval7 += valMask*roi8p[ baseRoi +7 ] ;
2236           }
2237
2238   // 8 pixel par thread --> global mem
2239   output[ idx   ] = outval0 ;
2240   output[ idx+1 ] = outval1 ;
2241   output[ idx+2 ] = outval2 ;
2242   output[ idx+3 ] = outval3 ;
2243   output[ idx+4 ] = outval4 ;
2244   output[ idx+5 ] = outval5 ;
2245   output[ idx+6 ] = outval6 ;
2246   output[ idx+7 ] = outval7 ;
2247 }
2248
2249
2250 __global__ void kernel_convoSepShx8pH(unsigned char *output, int i_dim, int j_dim, int r)
2251 {
2252   int idb, ic, jc, p;
2253   int L = 2*r+1 ;
2254   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 ;
2255   int bdimX = blockDim.x<<3 ;
2256   int tidX = threadIdx.x<<3 ;
2257
2258   // coordonnees absolues du point de base
2259   int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<3 ;
2260   int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
2261   int j0= __umul24(blockIdx.x,blockDim.x)<<3 ;
2262   int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image
2263
2264
2265   // chargement en smem
2266   int idrow = threadIdx.y*(bdimX+L-1) ;
2267
2268   extern __shared__ unsigned char roi8p[];
2269
2270   // bloc 0 (a gauche)
2271   for (p=0; p<8; p++)
2272         roi8p[  idrow  + tidX +p ] = tex2D(tex_img_inc, j-r+p  , i) ;
2273
2274   // a droite
2275   if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
2276         {
2277           roi8p[  idrow + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i ) ;
2278           roi8p[  idrow + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i ) ;
2279         }
2280
2281   __syncthreads();
2282
2283   // calculs de convolution
2284   // passe horizontale
2285   for (jc=0 ; jc<L ; jc++)
2286           {
2287                 int baseRoi = idrow + tidX +jc ;
2288                 float valMask = masque[ jc ] ;
2289                 outval0 += valMask*roi8p[  baseRoi    ] ;
2290                 outval1 += valMask*roi8p[  baseRoi +1 ] ;
2291                 outval2 += valMask*roi8p[  baseRoi +2 ] ;
2292                 outval3 += valMask*roi8p[  baseRoi +3 ] ;
2293                 outval4 += valMask*roi8p[  baseRoi +4 ] ;
2294                 outval5 += valMask*roi8p[  baseRoi +5 ] ;
2295                 outval6 += valMask*roi8p[  baseRoi +6 ] ;
2296                 outval7 += valMask*roi8p[  baseRoi +7 ] ;
2297           }
2298
2299   // 1 pixel par thread --> global mem
2300   output[ idx   ] = outval0 ;
2301   output[ idx+1 ] = outval1 ;
2302   output[ idx+2 ] = outval2 ;
2303   output[ idx+3 ] = outval3 ;
2304   output[ idx+4 ] = outval4 ;
2305   output[ idx+5 ] = outval5 ;
2306   output[ idx+6 ] = outval6 ;
2307   output[ idx+7 ] = outval7 ;
2308 }
2309
2310
2311 /*************************************************************************************************
2312  ***********************************************************************************************
2313
2314               FIN DES kERNELS de CONVOLUTION
2315
2316  ***********************************************************************************************
2317  *************************************************************************************************/
2318 // kernel de la libjacket
2319 // Exchange trick: Morgan McGuire, ShaderX 2008
2320 #define s2(a,b)            { float tmp = a; a = min(a,b); b = max(tmp,b); }
2321 #define mn3(a,b,c)         s2(a,b); s2(a,c);
2322 #define mx3(a,b,c)         s2(b,c); s2(a,c);
2323
2324 #define mnmx3(a,b,c)       mx3(a,b,c); s2(a,b);                               // 3 exchanges
2325 #define mnmx4(a,b,c,d)     s2(a,b); s2(c,d); s2(a,c); s2(b,d);                // 4 exchanges
2326 #define mnmx5(a,b,c,d,e)   s2(a,b); s2(c,d); mn3(a,c,e); mx3(b,d,e);          // 6 exchanges
2327 #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
2328
2329 #define IN(X,Y)  ((0 <= (X) && (X) < nx && 0 <= (Y) && (Y) < ny) ? d_in[(Y)*nx+(X)] : 0)
2330
2331 __global__ static void kernel_medjacket(int nx, int ny, unsigned short*d_out, unsigned short*d_in)
2332 {
2333   int x = __mul24(blockIdx.x , blockDim.x) + threadIdx.x;
2334   int y = __mul24(blockIdx.y , blockDim.y) + threadIdx.y;
2335
2336     // pull top six from image
2337     float v[6] = { IN(x-1, y-1), IN(x  , y-1), IN(x+1, y-1),
2338                    IN(x-1, y  ), IN(x  , y  ), IN(x+1, y  ) };
2339
2340     // with each pass, remove min and max values and add new value
2341     mnmx6(v[0], v[1], v[2], v[3], v[4], v[5]);
2342     v[5] = IN(x-1, y+1); // add new contestant
2343     mnmx5(v[1], v[2], v[3], v[4], v[5]);
2344     v[5] = IN(x  , y+1); // add new contestant
2345     mnmx4(v[2], v[3], v[4], v[5]);
2346     v[5] = IN(x+1, y+1); // add last contenstant
2347     mnmx3(v[3], v[4], v[5]);
2348
2349     // pick the middle one
2350     d_out[y*nx + x] = v[4];
2351 }
2352
2353
2354 /***************************************************************
2355  *   fonction de tri de 2 valeurs entieres (min en premier)
2356  ***************************************************************/
2357 __device__ inline void s(int* a, int* b)
2358 {
2359
2360   int tmp ;
2361   if (*a > *b)
2362         {
2363           tmp = *b ;
2364           *b = *a ;
2365           *a = tmp ;
2366           }
2367 }
2368
2369 __device__ inline void s(unsigned short * a, unsigned short* b)
2370 {
2371
2372   unsigned short tmp ;
2373   if (*a > *b)
2374         {
2375           tmp = *b ;
2376           *b = *a ;
2377           *a = tmp ;
2378           }
2379 }
2380
2381 __device__ inline void s(unsigned char * a, unsigned char * b)
2382 {
2383
2384   unsigned short tmp ;
2385   if (*a > *b)
2386         {
2387           tmp = *b ;
2388           *b = *a ;
2389           *a = tmp ;
2390           }
2391 }
2392
2393 /***************************************************************
2394  *   fonction de min-max d'un tableau de n valeurs entieres
2395  ***************************************************************/
2396 __device__ void minmaxN(unsigned short* v, int n)
2397 {
2398   for (int i=1; i< n; i++)
2399         s( v, v+i) ;
2400   for (int i=n-2; i>0; i--)
2401         s( v+i, v+n-1) ;
2402 }
2403
2404 /***************************************************************
2405  *   fonction de tri naif d'un tableau de n valeurs entieres
2406  ***************************************************************/
2407 __device__ void bubTriN(int * v, int n)
2408 {
2409   for (int i=0; i< n-1; i++)
2410         for (int j=i+1; j<n; j++)
2411           s( v+i, v+j) ;
2412 }
2413
2414
2415 /*****************************************************************************************
2416  *    MACROS pour tri min-max utilisant la fonction de tri s(a,b) ci-dessus
2417  *****************************************************************************************/
2418 #define min3(a, b, c) s(a, b); s(a, c);
2419 #define max3(a, b, c) s(b, c); s(a, c);
2420 #define minmax3(a, b, c) max3(a, b, c); s(a, b); // max 3
2421 #define minmax4(a, b, c, d) s(a, b); s(c, d); s(a, c); s(b, d); //max 4
2422 #define minmax5(a, b, c, d, e) s(a, b); s(c, d); min3(a, c, e); max3(b, d, e); //max 6
2423 #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
2424
2425 /***************************************************************
2426  *   fonction de tri de 9 valeurs entieres
2427  ***************************************************************/
2428 #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);
2429
2430
2431 __global__ void kernel_bubMedian3( unsigned short*output, int i_dim, int j_dim)
2432 {
2433
2434   // coordonnees absolues du point
2435   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
2436   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2437
2438   /**************************************************************************
2439    *             tri(s)
2440    **************************************************************************/
2441   int a0, a1, a2, a3, a4, a5, a6, a7, a8 ;
2442
2443   /********************************************************************************
2444    * affectation des valeurs
2445    ********************************************************************************/
2446   a0 = tex2D(tex_img_ins, j-1, i-1) ;
2447   a1 = tex2D(tex_img_ins, j  , i-1) ;
2448   a2 = tex2D(tex_img_ins, j+1, i-1) ;
2449   a3 = tex2D(tex_img_ins, j-1, i) ;
2450   a4 = tex2D(tex_img_ins, j  , i) ;
2451   a5 = tex2D(tex_img_ins, j+1, i) ;
2452   a6 = tex2D(tex_img_ins, j-1, i+1) ;
2453   a7 = tex2D(tex_img_ins, j  , i+1) ;
2454   a8 = tex2D(tex_img_ins, j+1, i+1) ;
2455
2456   //tri selection
2457   bubTriReg9(&a0, &a1, &a2, &a3, &a4, &a5, &a6, &a7, &a8);
2458
2459   //median au milieu !
2460   output[ __mul24(i, j_dim) +j ] = a4 ;
2461 }
2462
2463 __global__ void kernel_median3( unsigned short*output, int i_dim, int j_dim)
2464 {
2465
2466   // coordonnees absolues du point
2467   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
2468   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2469
2470   /**************************************************************************
2471    *             tri(s)
2472    **************************************************************************/
2473   int a0, a1, a2, a3, a4, a5 ;
2474
2475   /********************************************************************************
2476    * les six premieres valeurs (suffisant pour median 3x3 par forgetfull selection)
2477    ********************************************************************************/
2478   a0 = tex2D(tex_img_ins, j-1, i-1) ;
2479   a1 = tex2D(tex_img_ins, j, i-1) ;
2480   a2 = tex2D(tex_img_ins, j+1, i-1) ;
2481   a3 = tex2D(tex_img_ins, j-1, i) ;
2482   a4 = tex2D(tex_img_ins, j, i) ;
2483   a5 = tex2D(tex_img_ins, j+1, i) ;
2484
2485
2486   //min max aux extremites
2487   minmax6(&a0, &a1, &a2, &a3, &a4, &a5) ;
2488
2489   /********************************************
2490    * les deux valeurs suivantes aux extremites
2491    ********************************************/
2492   a5 = tex2D(tex_img_in, j-1, i+1) ;
2493
2494   minmax5(&a1, &a2, &a3, &a4, &a5) ;
2495
2496   /********************************************
2497    * la derniere valeur a la fin
2498    ********************************************/
2499
2500   a5 = tex2D(tex_img_ins, j, i+1) ;
2501
2502   minmax4(&a2, &a3, &a4, &a5) ;
2503
2504   a5 = tex2D(tex_img_ins, j+1, i+1) ;
2505   minmax3(&a3, &a4, &a5) ;
2506
2507
2508   //median au milieu !
2509   output[ __mul24(i, j_dim) +j ] = a4 ;
2510
2511 }
2512
2513 __global__ void kernel_median3Sh( unsigned short*output, int i_dim, int j_dim)
2514 {
2515   int idb, ic, jc ;
2516   int L = 3 ;
2517
2518   // coordonnees absolues du point de base
2519   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
2520   int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
2521   int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
2522
2523   // chargement en smem
2524   int idrow = threadIdx.y*(blockDim.x+2) ;
2525
2526   extern __shared__ int medRoi[];
2527
2528   // bloc 0 (en haut \E0 gauche)
2529   medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-1, i-1) ;
2530   // bloc 1 (en haut \E0 droite)...
2531   if ( threadIdx.x < L-1 ) //...ou plutot ce qu'il en manque
2532         medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-1, i-1 ) ;
2533   // bloc 2 ( en bas \E0 gauche)
2534   if ( threadIdx.y < L-1 )
2535         {
2536           idrow = (threadIdx.y+blockDim.y)*(blockDim.x+L-1) ;
2537           medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-1, i+blockDim.y-1 ) ;
2538           //bloc 4 ( en bas \E0 doite )...
2539           if ( threadIdx.x < L-1 ) //...ou ce qu'il en manque
2540                 medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-1, i+blockDim.y-1 ) ;
2541         }
2542   __syncthreads();
2543
2544   /**************************************************************************
2545    *             tri(s)
2546    **************************************************************************/
2547   int a0, a1, a2, a3, a4, a5 ;
2548
2549   /********************************************************************************
2550    * les six premieres valeurs (suffisant pour median 3x3 par forgetfull selection)
2551    ********************************************************************************/
2552   // l'index dans medRoi est du type (threadIdx.y+ic)*(blockDim.X+L-1)+ (threadIdx.x+jc)
2553   // ou ic,jc sont les coordonnees du point dans le noyau 3x3
2554   a0 = medRoi[ (threadIdx.y)*(blockDim.x+L-1)+ (threadIdx.x+1) ] ;
2555   a1 = medRoi[ (threadIdx.y)*(blockDim.x+L-1)+ (threadIdx.x+2)  ] ;
2556   a2 = medRoi[ (threadIdx.y+1)*(blockDim.x+L-1)+ (threadIdx.x+1)] ;
2557   a3 = medRoi[ (threadIdx.y+1)*(blockDim.x+L-1)+ (threadIdx.x+2)] ;
2558   a4 = medRoi[ (threadIdx.y+2)*(blockDim.x+L-1)+ (threadIdx.x+1)] ;
2559   a5 = medRoi[ (threadIdx.y+2)*(blockDim.x+L-1)+ (threadIdx.x+2)] ;
2560
2561
2562   //min max aux extremites
2563   minmax6(&a0, &a1, &a2, &a3, &a4, &a5) ;
2564
2565   /********************************************
2566    * les deux valeurs suivantes aux extremites
2567    ********************************************/
2568   a5 = medRoi[ (threadIdx.y)*(blockDim.x+L-1)+ (threadIdx.x)] ;
2569
2570   minmax5(&a1, &a2, &a3, &a4, &a5) ;
2571
2572   /********************************************
2573    * la derniere valeur a la fin
2574    ********************************************/
2575
2576   a5 = medRoi[ (threadIdx.y+1)*(blockDim.x+L-1)+ (threadIdx.x)] ;
2577
2578   minmax4(&a2, &a3, &a4, &a5) ;
2579
2580   a5 = medRoi[ (threadIdx.y+2)*(blockDim.x+L-1)+ (threadIdx.x)] ;
2581   minmax3(&a3, &a4, &a5) ;
2582
2583
2584   //median au milieu !
2585   output[ __mul24(i, j_dim) +j ] = a4 ;
2586
2587 }
2588
2589 /************************************************
2590  * version qui gere 2 pxels par thread
2591  * Objectif plus déebit en cachant la latence
2592  * minimiser les acces en texture
2593  * defaut = plus de registres
2594  ************************************************/
2595 __global__ void kernel_median3_2pix( unsigned short*output, int i_dim, int j_dim)
2596 {
2597
2598   // coordonnees absolues du point
2599   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
2600   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2601
2602   /**************************************************************************
2603    *             tri(s)
2604    **************************************************************************/
2605   int a0, a1, a2, a3, a4, a5 ;
2606   int b0, b1, b2, b3, b4, b5 ;
2607
2608   /********************************************************************************
2609    * les six premieres valeurs (suffisant pour median 3x3 par forgetfull selection)
2610    ********************************************************************************/
2611   a0 = tex2D(tex_img_ins, j  , i-1) ;
2612   a1 = tex2D(tex_img_ins, j+1, i-1) ;
2613   a2 = tex2D(tex_img_ins, j  , i  ) ;
2614   a3 = tex2D(tex_img_ins, j+1, i  ) ;
2615   a4 = tex2D(tex_img_ins, j  , i+1) ;
2616   a5 = tex2D(tex_img_ins, j+1, i+1) ;
2617
2618
2619   //min max aux extremites
2620   minmax6(&a0, &a1, &a2, &a3, &a4, &a5) ;
2621
2622   b0=a0; b1=a1; b2=a2; b3=a3; b4=a4; b5=a5;
2623
2624   /********************************************
2625    * la valeur suivante au bout
2626    ********************************************/
2627   a5 = tex2D(tex_img_ins, j-1, i  ) ;
2628   b5 = tex2D(tex_img_ins, j+2, i  ) ;
2629
2630   minmax5(&a1, &a2, &a3, &a4, &a5) ;
2631   minmax5(&b1, &b2, &b3, &b4, &b5) ;
2632
2633   /********************************************
2634    * la derniere valeur a la fin
2635    ********************************************/
2636   a5 = tex2D(tex_img_ins, j-1, i-1) ;
2637   b5 = tex2D(tex_img_ins, j+2, i-1) ;
2638
2639   minmax4(&a2, &a3, &a4, &a5) ;
2640   minmax4(&b2, &b3, &b4, &b5) ;
2641
2642   a5 = tex2D(tex_img_ins, j-1, i+1) ;
2643   b5 = tex2D(tex_img_ins, j+2, i+1) ;
2644
2645   minmax3(&a3, &a4, &a5) ;
2646   minmax3(&b3, &b4, &b5) ;
2647
2648   //median au milieu !
2649   output[ __mul24(i, j_dim) +j ] = a4 ;
2650   output[ __mul24(i, j_dim) +j+1 ] = b4 ;
2651
2652 }
2653
2654 __global__ void kernel_median3_2pix8( unsigned char  *output, int i_dim, int j_dim)
2655 {
2656
2657   // coordonnees absolues du point
2658   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
2659   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2660
2661   /**************************************************************************
2662    *             tri(s)
2663    **************************************************************************/
2664   int a0, a1, a2, a3, a4, a5 ;
2665   int b0, b1, b2, b3, b4, b5 ;
2666
2667   /********************************************************************************
2668    * les six premieres valeurs (suffisant pour median 3x3 par forgetfull selection)
2669    ********************************************************************************/
2670   a0 = tex2D(tex_img_inc, j  , i-1) ;
2671   a1 = tex2D(tex_img_inc, j+1, i-1) ;
2672   a2 = tex2D(tex_img_inc, j  , i  ) ;
2673   a3 = tex2D(tex_img_inc, j+1, i  ) ;
2674   a4 = tex2D(tex_img_inc, j  , i+1) ;
2675   a5 = tex2D(tex_img_inc, j+1, i+1) ;
2676
2677
2678   //min max aux extremites
2679   minmax6(&a0, &a1, &a2, &a3, &a4, &a5) ;
2680
2681   b0=a0; b1=a1; b2=a2; b3=a3; b4=a4; b5=a5;
2682
2683   /********************************************
2684    * la valeur suivante au bout
2685    ********************************************/
2686   a5 = tex2D(tex_img_inc, j-1, i  ) ;
2687   b5 = tex2D(tex_img_inc, j+2, i  ) ;
2688
2689   minmax5(&a1, &a2, &a3, &a4, &a5) ;
2690   minmax5(&b1, &b2, &b3, &b4, &b5) ;
2691
2692   /********************************************
2693    * la derniere valeur a la fin
2694    ********************************************/
2695   a5 = tex2D(tex_img_inc, j-1, i-1) ;
2696   b5 = tex2D(tex_img_inc, j+2, i-1) ;
2697
2698   minmax4(&a2, &a3, &a4, &a5) ;
2699   minmax4(&b2, &b3, &b4, &b5) ;
2700
2701   a5 = tex2D(tex_img_inc, j-1, i+1) ;
2702   b5 = tex2D(tex_img_inc, j+2, i+1) ;
2703
2704   minmax3(&a3, &a4, &a5) ;
2705   minmax3(&b3, &b4, &b5) ;
2706
2707   //median au milieu !
2708   output[ __mul24(i, j_dim) +j ] = a4 ;
2709   output[ __mul24(i, j_dim) +j+1 ] = b4 ;
2710
2711 }
2712
2713
2714 #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
2715 #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
2716 #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
2717 #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
2718 #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
2719 #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
2720 #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 !!
2721 #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
2722 #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
2723 #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
2724
2725 #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
2726 #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
2727 #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
2728 #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
2729 #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 !!!
2730 #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
2731 #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
2732 #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
2733 #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
2734 #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
2735 #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
2736 #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
2737 #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
2738
2739
2740 __global__ void kernel_median5( unsigned short*output, int i_dim, int j_dim)
2741 {
2742
2743   // coordonnees absolues du point
2744   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
2745   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2746
2747   /**************************************************************************
2748    *             tri(s)
2749    **************************************************************************/
2750   int a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13 ;
2751
2752   /********************************************************************************
2753    * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
2754    ********************************************************************************/
2755   //premiere ligne
2756   a0  = tex2D(tex_img_ins, j-2, i-2) ;
2757   a1  = tex2D(tex_img_ins, j-1, i-2) ;
2758   a2  = tex2D(tex_img_ins, j  , i-2) ;
2759   a3  = tex2D(tex_img_ins, j+1, i-2) ;
2760   a4  = tex2D(tex_img_ins, j+2, i-2) ;
2761   //deuxieme ligne
2762   a5  = tex2D(tex_img_ins, j-2, i-1) ;
2763   a6  = tex2D(tex_img_ins, j-1, i-1) ;
2764   a7  = tex2D(tex_img_ins, j  , i-1) ;
2765   a8  = tex2D(tex_img_ins, j+1, i-1) ;
2766   a9  = tex2D(tex_img_ins, j+2, i-1) ;
2767   //les 4 premiers de la troisieme ligne
2768   a10 = tex2D(tex_img_ins, j-2, i) ;
2769   a11 = tex2D(tex_img_ins, j-1, i) ;
2770   a12 = tex2D(tex_img_ins, j  , i) ;
2771   a13 = tex2D(tex_img_ins, j+1, i) ;
2772
2773   //min max aux extremites
2774   minmax14(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2775
2776   //chargement valeurs suivante (15)
2777   a13 = tex2D(tex_img_ins, j+2, i);
2778   //minmax aux extremites
2779   minmax13(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2780
2781
2782   //chargement valeur suivante (16)
2783   a13 = tex2D(tex_img_ins, j-2, i+1);
2784   //minmax aux extremites
2785   minmax12(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2786
2787   //chargement valeur suivante (17)
2788   a13 = tex2D(tex_img_ins, j-1, i+1);
2789   //minmax aux extremites
2790   minmax11(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2791
2792   //chargement valeur suivante (18)
2793   a13 = tex2D(tex_img_ins, j  , i+1);
2794   //minmax aux extremites
2795   minmax10(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2796
2797   //chargement valeur suivante (19)
2798   a13 = tex2D(tex_img_ins, j+1, i+1);
2799   //minmax aux extremites
2800   minmax9(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2801
2802   //chargement valeur suivante (20)
2803   a13 = tex2D(tex_img_ins, j+2, i+1);
2804   //minmax aux extremites
2805   minmax8(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2806
2807   //chargement valeur suivante (21)
2808   a13 = tex2D(tex_img_ins, j-2, i+2);
2809   //minmax aux extremites
2810   minmax7(&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2811
2812   //chargement valeur suivante (22)
2813   a13 = tex2D(tex_img_ins, j-1, i+2);
2814   //minmax aux extremites
2815   minmax6(&a8,&a9,&a10,&a11,&a12,&a13);
2816
2817   //chargement valeur suivante (23)
2818   a13 = tex2D(tex_img_ins, j  , i+2);
2819   //minmax aux extremites
2820   minmax5(&a9,&a10,&a11,&a12,&a13);
2821
2822   //chargement valeur suivante (24)
2823   a13 = tex2D(tex_img_ins, j+1  , i+2);
2824   //minmax aux extremites
2825   minmax4(&a10,&a11,&a12,&a13);
2826
2827   //chargement valeur suivante (25)
2828   a13 = tex2D(tex_img_ins, j+2  , i+2);
2829   //minmax aux extremites
2830   minmax3(&a11,&a12,&a13);
2831
2832
2833   //median au milieu !
2834   output[ __mul24(i, j_dim) +j ] = a12 ;
2835
2836 }
2837
2838 __global__ void kernel_median5_sh( unsigned short*output, int i_dim, int j_dim)
2839 {
2840
2841   // coordonnees absolues du point
2842   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
2843   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2844   int L = 5 ;
2845   int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
2846
2847   // chargement en smem
2848   int idrow = threadIdx.y*(blockDim.x+4) ;
2849
2850   extern __shared__ int medRoi[];
2851
2852   // bloc 0 (en haut \E0 gauche)
2853   medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-2, i-2) ;
2854   // bloc 1 (en haut \E0 droite)...
2855   if ( threadIdx.x < L-1 ) //...ou plutot ce qu'il en manque
2856         medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-2, i-2 ) ;
2857   // bloc 2 ( en bas \E0 gauche)
2858   if ( threadIdx.y < L-1 )
2859         {
2860           idrow = (threadIdx.y+blockDim.y)*(blockDim.x+L-1) ;
2861           medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-2, i+blockDim.y-2 ) ;
2862           //bloc 4 ( en bas \E0 doite )...
2863           if ( threadIdx.x < L-1 ) //...ou ce qu'il en manque
2864                 medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-2, i+blockDim.y-2 ) ;
2865         }
2866   __syncthreads();
2867
2868   /**************************************************************************
2869    *             tri(s)
2870    **************************************************************************/
2871   int a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13 ;
2872
2873   /********************************************************************************
2874    * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
2875    ********************************************************************************/
2876   // l'index dans medRoi est du type (threadIdx.y+ic)*(blockDim.x+L-1)+ (threadIdx.x+jc)
2877   // ou ic,jc sont les coordonnees du point dans le noya
2878   //premiere ligne
2879   a0  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i-2) ;
2880   a1  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i-2) ;
2881   a2  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i-2) ;
2882   a3  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i-2) ;
2883   a4  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i-2) ;
2884   //deuxieme ligne
2885   a5  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i-1) ;
2886   a6  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i-1) ;
2887   a7  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i-1) ;
2888   a8  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i-1) ;
2889   a9  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i-1) ;
2890   //les 4 premiers de la troisieme ligne
2891   a10 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i) ;
2892   a11 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i) ;
2893   a12 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i) ;
2894   a13 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i) ;
2895
2896   //min max aux extremites
2897   minmax14(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2898
2899   //chargement valeurs suivante (15)
2900   a13 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i);
2901   //minmax aux extremites
2902   minmax13(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2903
2904
2905   //chargement valeur suivante (16)
2906   a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i+1);
2907   //minmax aux extremites
2908   minmax12(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2909
2910   //chargement valeur suivante (17)
2911   a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i+1);
2912   //minmax aux extremites
2913   minmax11(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2914
2915   //chargement valeur suivante (18)
2916   a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i+1);
2917   //minmax aux extremites
2918   minmax10(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2919
2920   //chargement valeur suivante (19)
2921   a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i+1);
2922   //minmax aux extremites
2923   minmax9(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2924
2925   //chargement valeur suivante (20)
2926   a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i+1);
2927   //minmax aux extremites
2928   minmax8(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2929
2930   //chargement valeur suivante (21)
2931   a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i+2);
2932   //minmax aux extremites
2933   minmax7(&a7,&a8,&a9,&a10,&a11,&a12,&a13);
2934
2935   //chargement valeur suivante (22)
2936   a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i+2);
2937   //minmax aux extremites
2938   minmax6(&a8,&a9,&a10,&a11,&a12,&a13);
2939
2940   //chargement valeur suivante (23)
2941   a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i+2);
2942   //minmax aux extremites
2943   minmax5(&a9,&a10,&a11,&a12,&a13);
2944
2945   //chargement valeur suivante (24)
2946   a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1  , i+2);
2947   //minmax aux extremites
2948   minmax4(&a10,&a11,&a12,&a13);
2949
2950   //chargement valeur suivante (25)
2951   a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2  , i+2);
2952   //minmax aux extremites
2953   minmax3(&a11,&a12,&a13);
2954
2955
2956   //median au milieu !
2957   output[ __mul24(i, j_dim) +j ] = a12 ;
2958
2959 }
2960
2961
2962 __global__ void kernel_median5_sh_loc( unsigned short*output, int i_dim, int j_dim)
2963 {
2964
2965   // coordonnees absolues du point
2966   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
2967   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
2968   // int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
2969
2970   // chargement en smem
2971   int idrow = threadIdx.y*(blockDim.x+4) ;
2972
2973   extern __shared__ int medRoi[];
2974
2975   // bloc 0 (en haut \E0 gauche)
2976   medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-2, i-2) ;
2977   // bloc 1 (en haut \E0 droite)...
2978   if ( threadIdx.x < 4 ) //...ou plutot ce qu'il en manque
2979         medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-2, i-2 ) ;
2980   // bloc 2 ( en bas \E0 gauche)
2981   if ( threadIdx.y < 4 )
2982         {
2983           idrow = (threadIdx.y+blockDim.y)*(blockDim.x+4) ;
2984           medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-2, i+blockDim.y-2 ) ;
2985           //bloc 4 ( en bas \E0 doite )...
2986           if ( threadIdx.x < 4 ) //...ou ce qu'il en manque
2987                 medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-2, i+blockDim.y-2 ) ;
2988         }
2989   __syncthreads();
2990
2991   /**************************************************************************
2992    *             tri(s)
2993    **************************************************************************/
2994   /********************************************************************************
2995    * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
2996    ********************************************************************************/
2997   // l'index dans medRoi est du type (threadIdx.y+ic)*(blockDim.x+L-1)+ (threadIdx.x+jc)
2998   // ou ic,jc sont les coordonnees du point dans le noyau
2999
3000   int a[] = { medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x ) ],
3001                           medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x +1) ],
3002                           medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x +2) ],
3003                           medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x +3) ],
3004                           medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x +4) ],
3005                           medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x ) ],
3006                           medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x +1) ],
3007                           medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x +2) ],
3008                           medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x +3) ],
3009                           medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x +4) ],
3010                           medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x ) ],
3011                           medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x +1) ],
3012                           medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x +2) ],
3013                           medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x +3) ]};
3014
3015
3016   //min max aux extremites
3017   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]);
3018
3019   //chargement valeurs suivante (15)
3020   a[13] = medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i);
3021   //minmax aux extremites
3022   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]);
3023
3024
3025   //chargement valeur suivante (16)
3026   a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i+1);
3027   //minmax aux extremites
3028   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]);
3029
3030   //chargement valeur suivante (17)
3031   a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i+1);
3032   //minmax aux extremites
3033   minmax11(&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
3034
3035   //chargement valeur suivante (18)
3036   a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i+1);
3037   //minmax aux extremites
3038   minmax10(&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
3039
3040   //chargement valeur suivante (19)
3041   a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i+1);
3042   //minmax aux extremites
3043   minmax9(&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
3044
3045   //chargement valeur suivante (20)
3046   a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i+1);
3047   //minmax aux extremites
3048   minmax8(&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
3049
3050   //chargement valeur suivante (21)
3051   a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i+2);
3052   //minmax aux extremites
3053   minmax7(&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
3054
3055   //chargement valeur suivante (22)
3056   a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i+2);
3057   //minmax aux extremites
3058   minmax6(&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
3059
3060   //chargement valeur suivante (23)
3061   a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i+2);
3062   //minmax aux extremites
3063   minmax5(&a[9],&a[10],&a[11],&a[12],&a[13]);
3064
3065   //chargement valeur suivante (24)
3066   a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1  , i+2);
3067   //minmax aux extremites
3068   minmax4(&a[10],&a[11],&a[12],&a[13]);
3069
3070   //chargement valeur suivante (25)
3071   a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2  , i+2);
3072   //minmax aux extremites
3073   minmax3(&a[11],&a[12],&a[13]);
3074
3075
3076   //median au milieu !
3077   output[ __mul24(i, j_dim) +j ] = a[12] ;
3078
3079 }
3080
3081 __global__ void kernel_median5_sh_loc_2p( unsigned short*output, int i_dim, int j_dim)
3082 {
3083
3084   // coordonnees absolues du point
3085   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
3086   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
3087   // int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
3088   int bdimX = blockDim.x<<1 ;
3089   int tidX = threadIdx.x<<1 ;
3090
3091   // chargement en smem
3092   int idrow = threadIdx.y*(bdimX+2*2) ;
3093
3094   extern __shared__ int medRoi[];
3095
3096   // bloc 0 (en haut \E0 gauche)
3097   medRoi[ idrow  + tidX ] = tex2D(tex_img_ins, j-2, i-2) ;
3098   medRoi[ idrow  + tidX +1] = tex2D(tex_img_ins, j-1, i-2) ;
3099   // bloc 1 (en haut \E0 droite)...
3100   if ( tidX < 4 ) //...ou plutot ce qu'il en manque
3101         {
3102           medRoi[ idrow + tidX + bdimX ] = tex2D( tex_img_ins, j+bdimX-2, i-2 ) ;
3103           medRoi[ idrow + tidX + bdimX +1] = tex2D( tex_img_ins, j+bdimX-1, i-2 ) ;
3104         }
3105   // bloc 2 ( en bas \E0 gauche)
3106   if ( threadIdx.y < 4 )
3107         {
3108           idrow = (threadIdx.y+blockDim.y)*(bdimX+4) ;
3109           medRoi[ idrow + tidX ] = tex2D( tex_img_ins, j-2, i+blockDim.y-2 ) ;
3110           medRoi[ idrow + tidX +1] = tex2D( tex_img_ins, j-1, i+blockDim.y-2 ) ;
3111           //bloc 4 ( en bas \E0 doite )...
3112           if ( tidX < 4 ) //...ou ce qu'il en manque
3113                 {
3114                   medRoi[ idrow + tidX + bdimX ] = tex2D( tex_img_ins, j+bdimX-2, i+blockDim.y-2 ) ;
3115                   medRoi[ idrow + tidX + bdimX +1] = tex2D( tex_img_ins, j+bdimX-1, i+blockDim.y-2 ) ;
3116                 }
3117         }
3118   __syncthreads();
3119
3120   /**************************************************************************
3121    *             tri(s)
3122    **************************************************************************/
3123   /********************************************************************************
3124    * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
3125    ********************************************************************************/
3126   // l'index dans medRoi est du type (threadIdx.y+ic)*(bdimX+L-1)+ (tidX+jc)
3127   // ou ic,jc sont les coordonnees du point dans le noyau
3128
3129   //on commence le tri avec les pixels du voisinage commun aux 2 pixels \E0 calculer
3130
3131   int a[] = { medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +1) ],
3132                           medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +2) ],
3133                           medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +3) ],
3134                           medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +4) ],
3135
3136                           medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +1) ],
3137                           medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +2) ],
3138                           medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +3) ],
3139                           medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +4) ],
3140
3141                           medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +1) ],
3142                           medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +2) ],
3143                           medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +3) ],
3144                           medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +4) ],
3145
3146                           medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +1) ],
3147                           medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +2) ]
3148   };
3149
3150   //min max aux extremites
3151   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]);
3152
3153   //chargement valeurs suivante (15)
3154   a[13] = medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +3) ] ;
3155   //minmax aux extremites
3156   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]);
3157
3158
3159   //chargement valeur suivante (16)
3160   a[13] = medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +4) ] ;
3161   //minmax aux extremites
3162   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]);
3163
3164   //chargement valeur suivante (17)
3165   a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +1) ] ;
3166   //minmax aux extremites
3167   minmax11(&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
3168
3169   //chargement valeur suivante (18)
3170   a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +2) ] ;
3171   //minmax aux extremites
3172   minmax10(&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
3173
3174   //chargement valeur suivante (19)
3175   a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +3) ] ;
3176   //minmax aux extremites
3177   minmax9(&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
3178
3179   //chargement valeur suivante (20)
3180   a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +4) ] ;
3181   //minmax aux extremites
3182   minmax8(&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
3183
3184   // fin des pixels communs. la suite doit se faire en 2 tris distincts.
3185   int b[] = { a[7], a[8], a[9], a[10], a[11], a[12], a[13]};
3186
3187   //chargement valeur suivante (21)
3188   a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX ) ] ;
3189   b[6 ] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +5) ] ;
3190   //minmax aux extremites
3191   minmax7(&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
3192   minmax7(&b[0],&b[1],&b[2],&b[3],&b[4],&b[5],&b[6]);
3193
3194   //chargement valeur suivante (22)
3195   a[13] = medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX ) ] ;
3196   b[6 ] = medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +5) ] ;
3197   //minmax aux extremites
3198   minmax6(&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
3199   minmax6(&b[1],&b[2],&b[3],&b[4],&b[5],&b[6]);
3200
3201   //chargement valeur suivante (23)
3202   a[13] = medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX ) ] ;
3203   b[6 ] = medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +5) ] ;
3204   //minmax aux extremites
3205   minmax5(&a[9],&a[10],&a[11],&a[12],&a[13]);
3206   minmax5(&b[2],&b[3],&b[4],&b[5],&b[6]);
3207
3208   //chargement valeur suivante (24)
3209   a[13] = medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX ) ] ;
3210   b[6 ] = medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +5) ] ;
3211
3212   //minmax aux extremites
3213   minmax4(&a[10],&a[11],&a[12],&a[13]);
3214   minmax4(&b[3],&b[4],&b[5],&b[6]);
3215
3216   //chargement valeur suivante (25)
3217   a[13] = medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX ) ] ;
3218   b[6 ] = medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +5) ] ;
3219   //minmax aux extremites
3220   minmax3(&a[11],&a[12],&a[13]);
3221   minmax3(&b[4],&b[5],&b[6]);
3222
3223   //median au milieu !
3224   output[ __mul24(i, j_dim) +j ] = a[12] ;
3225   output[ __mul24(i, j_dim) +j+1 ] = b[5] ;
3226
3227 }
3228
3229
3230 __global__ void kernel_medianR( unsigned short * output, int i_dim, int j_dim, int r)
3231 {
3232   // coordonnees absolues du point
3233   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
3234   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
3235   unsigned short cpt ;
3236   unsigned short ic, jc ;
3237
3238   unsigned short hist[256] ;
3239   for (ic =0; ic<256; ic++) hist[ic]=0; // init histogramme
3240   //for (ic =0; ic<25; ic++) hist[ic]=1; // init histogramme
3241
3242   //generation histogramme
3243   for(ic=i-r; ic<=i+r; ic++ )
3244         for(jc=j-r; jc<=j+r; jc++)
3245           hist[tex2D(tex_img_ins, jc, ic)]++ ;
3246
3247   //parcours histogramme
3248   cpt = 0 ;
3249   for(ic=0; ic<256; ic++)
3250         {
3251           cpt += hist[ic] ;
3252           //selection de la valeur du percentile (ici 50%=SUM/2)
3253           if ( cpt > ((2*r+1)*(2*r+1))>>1 ) break ;
3254         }
3255
3256   output[ __mul24(i, j_dim) +j ] =  ic ;
3257 }
3258
3259 __global__ void kernel_medianRSH( unsigned short * output, int i_dim, int j_dim, int r)
3260 {
3261   // coordonnees absolues du point
3262   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
3263   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
3264   unsigned short cpt ;
3265   unsigned short ic, jc ;
3266
3267   // chargement en smem
3268   int idrow = threadIdx.y*(blockDim.x+2*r) ;
3269
3270   extern __shared__ int medRoi[];
3271
3272   // bloc 0 (en haut \E0 gauche)
3273   medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-r, i-r) ;
3274   // bloc 1 (en haut \E0 droite)...
3275   if ( threadIdx.x < 2*r ) //...ou plutot ce qu'il en manque
3276         medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-r, i-r ) ;
3277   // bloc 2 ( en bas \E0 gauche)
3278   if ( threadIdx.y < 2*r )
3279         {
3280           idrow = (threadIdx.y+blockDim.y)*(blockDim.x+2*r) ;
3281           medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-r, i+blockDim.y-r ) ;
3282           //bloc 4 ( en bas \E0 doite )...
3283           if ( threadIdx.x < 2*r ) //...ou ce qu'il en manque
3284                 medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-r, i+blockDim.y-r ) ;
3285         }
3286   __syncthreads();
3287
3288   int hist[256] ;
3289   for (ic =0; ic<256; ic++) hist[ic]=0; // init histogramme
3290
3291   //generation histogramme
3292   for(ic=0; ic<=2*r; ic++ )
3293         for(jc=0; jc<=2*r; jc++)
3294           hist[medRoi[(threadIdx.y+ic)*(blockDim.x+2*r)+ (threadIdx.x+jc)]]++ ;
3295
3296   //parcours histogramme
3297   cpt = 0 ;
3298   for(ic=0; ic<256; ic++)
3299         {
3300           cpt += hist[ic] ;
3301           //selection de la valeur du percentile (ici 50%=SUM/2)
3302           if ( cpt > (((2*r+1)*(2*r+1))>>1) ) break ;
3303         }
3304
3305   output[ __mul24(i, j_dim) +j ] =  ic ;
3306 }
3307
3308
3309 __global__ void kernel_median5_2pix( unsigned short*output, int i_dim, int j_dim)
3310 {
3311
3312   // coordonnees absolues du point
3313   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
3314   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
3315
3316   /**************************************************************************
3317    *             tri(s)
3318    **************************************************************************/
3319   int a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13 ;
3320   int b7, b8, b9, b10, b11, b12, b13 ;
3321
3322   /********************************************************************************
3323    * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
3324    ********************************************************************************/
3325   //premiere ligne
3326   a0  = tex2D(tex_img_ins, j-1, i-2) ;
3327   a1  = tex2D(tex_img_ins, j  , i-2) ;
3328   a2  = tex2D(tex_img_ins, j+1, i-2) ;
3329   a3  = tex2D(tex_img_ins, j+2, i-2) ;
3330   //deuxieme ligne
3331   a4  = tex2D(tex_img_ins, j-1, i-1) ;
3332   a5  = tex2D(tex_img_ins, j  , i-1) ;
3333   a6  = tex2D(tex_img_ins, j+1, i-1) ;
3334   a7  = tex2D(tex_img_ins, j+2, i-1) ;
3335   //troisieme ligne
3336   a8  = tex2D(tex_img_ins, j-1, i) ;
3337   a9  = tex2D(tex_img_ins, j  , i) ;
3338   a10 = tex2D(tex_img_ins, j+1, i) ;
3339   a11 = tex2D(tex_img_ins, j+2, i) ;
3340   //les 2 premiers de la 4eme ligne
3341   a12 = tex2D(tex_img_ins, j-1, i+1) ;
3342   a13 = tex2D(tex_img_ins, j  , i+1) ;
3343
3344   //min max aux extremites
3345   minmax14(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3346
3347   //chargement valeurs suivante (15)
3348   a13 = tex2D(tex_img_ins, j+1, i+1);
3349   //minmax aux extremites
3350   minmax13(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3351
3352   //chargement valeur suivante (16)
3353   a13 = tex2D(tex_img_ins, j+2, i+1);
3354   //minmax aux extremites
3355   minmax12(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3356
3357   //chargement valeur suivante (17)
3358   a13 = tex2D(tex_img_ins, j-1, i+2);
3359   //minmax aux extremites
3360   minmax11(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3361
3362   //chargement valeur suivante (18)
3363   a13 = tex2D(tex_img_ins, j  , i+2);
3364   //minmax aux extremites
3365   minmax10(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3366
3367   //chargement valeur suivante (19)
3368   a13 = tex2D(tex_img_ins, j+1, i+2);
3369   //minmax aux extremites
3370   minmax9(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3371
3372   //chargement valeur suivante (20)
3373   a13 = tex2D(tex_img_ins, j+2, i+2);
3374   //minmax aux extmites
3375   minmax8(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3376
3377   // fin des pixels voisins communs deux pixels centraux
3378   b7=a7; b8=a8; b9=a9; b10=a10; b11=a11; b12=a12; b13=a13;
3379
3380   //chargement valeur suivante (21)
3381   a13 = tex2D(tex_img_ins, j-2, i-2);
3382   b13 = tex2D(tex_img_ins, j+3, i-2);
3383   //minmax aux extremites
3384   minmax7(&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3385   minmax7(&b7,&b8,&b9,&b10,&b11,&b12,&b13);
3386
3387   //chargement valeur suivante (22)
3388   a13 = tex2D(tex_img_ins, j-2, i-1);
3389   b13 = tex2D(tex_img_ins, j+3, i-1);
3390   //minmax aux extremites
3391   minmax6(&a8,&a9,&a10,&a11,&a12,&a13);
3392   minmax6(&b8,&b9,&b10,&b11,&b12,&b13);
3393
3394   //chargement valeur suivante (23)
3395   a13 = tex2D(tex_img_ins, j-2, i  );
3396   b13 = tex2D(tex_img_ins, j+3, i  );
3397   //minmax aux extremites
3398   minmax5(&a9,&a10,&a11,&a12,&a13);
3399   minmax5(&b9,&b10,&b11,&b12,&b13);
3400
3401   //chargement valeur suivante (24)
3402   a13 = tex2D(tex_img_ins, j-2, i+1);
3403   b13 = tex2D(tex_img_ins, j+3, i+1);
3404   //minmax aux extremites
3405   minmax4(&a10,&a11,&a12,&a13);
3406   minmax4(&b10,&b11,&b12,&b13);
3407
3408   //chargement valeur suivante (25)
3409   a13 = tex2D(tex_img_ins, j-2, i+2);
3410   b13 = tex2D(tex_img_ins, j+3, i+2);
3411   //minmax aux extremites
3412   minmax3(&a11,&a12,&a13);
3413   minmax3(&b11,&b12,&b13);
3414
3415   //median au milieu !
3416   output[ __mul24(i, j_dim) +j ] = a12 ;
3417   output[ __mul24(i, j_dim) +j+1 ] = b12 ;
3418
3419 }
3420
3421
3422 __global__ void kernel_median5_2pix( unsigned char  *output, int i_dim, int j_dim)
3423 {
3424
3425   // coordonnees absolues du point
3426   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
3427   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
3428
3429   /**************************************************************************
3430    *             tri(s)
3431    **************************************************************************/
3432   int a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13 ;
3433   int b7, b8, b9, b10, b11, b12, b13 ;
3434
3435   /********************************************************************************
3436    * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
3437    ********************************************************************************/
3438   //premiere ligne
3439   a0  = tex2D(tex_img_inc, j-1, i-2) ;
3440   a1  = tex2D(tex_img_inc, j  , i-2) ;
3441   a2  = tex2D(tex_img_inc, j+1, i-2) ;
3442   a3  = tex2D(tex_img_inc, j+2, i-2) ;
3443   //deuxieme ligne
3444   a4  = tex2D(tex_img_inc, j-1, i-1) ;
3445   a5  = tex2D(tex_img_inc, j  , i-1) ;
3446   a6  = tex2D(tex_img_inc, j+1, i-1) ;
3447   a7  = tex2D(tex_img_inc, j+2, i-1) ;
3448   //troisieme ligne
3449   a8  = tex2D(tex_img_inc, j-1, i) ;
3450   a9  = tex2D(tex_img_inc, j  , i) ;
3451   a10 = tex2D(tex_img_inc, j+1, i) ;
3452   a11 = tex2D(tex_img_inc, j+2, i) ;
3453   //les 2 premiers de la 4eme ligne
3454   a12 = tex2D(tex_img_inc, j-1, i+1) ;
3455   a13 = tex2D(tex_img_inc, j  , i+1) ;
3456
3457   //min max aux extremites
3458   minmax14(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3459
3460   //chargement valeurs suivante (15)
3461   a13 = tex2D(tex_img_inc, j+1, i+1);
3462   //minmax aux extremites
3463   minmax13(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3464
3465   //chargement valeur suivante (16)
3466   a13 = tex2D(tex_img_inc, j+2, i+1);
3467   //minmax aux extremites
3468   minmax12(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3469
3470   //chargement valeur suivante (17)
3471   a13 = tex2D(tex_img_inc, j-1, i+2);
3472   //minmax aux extremites
3473   minmax11(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3474
3475   //chargement valeur suivante (18)
3476   a13 = tex2D(tex_img_inc, j  , i+2);
3477   //minmax aux extremites
3478   minmax10(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3479
3480   //chargement valeur suivante (19)
3481   a13 = tex2D(tex_img_inc, j+1, i+2);
3482   //minmax aux extremites
3483   minmax9(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3484
3485   //chargement valeur suivante (20)
3486   a13 = tex2D(tex_img_inc, j+2, i+2);
3487   //minmax aux extmites
3488   minmax8(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3489
3490   // fin des pixels voisins communs deux pixels centraux
3491   b7=a7; b8=a8; b9=a9; b10=a10; b11=a11; b12=a12; b13=a13;
3492
3493   //chargement valeur suivante (21)
3494   a13 = tex2D(tex_img_inc, j-2, i-2);
3495   b13 = tex2D(tex_img_inc, j+3, i-2);
3496   //minmax aux extremites
3497   minmax7(&a7,&a8,&a9,&a10,&a11,&a12,&a13);
3498   minmax7(&b7,&b8,&b9,&b10,&b11,&b12,&b13);
3499
3500   //chargement valeur suivante (22)
3501   a13 = tex2D(tex_img_inc, j-2, i-1);
3502   b13 = tex2D(tex_img_inc, j+3, i-1);
3503   //minmax aux extremites
3504   minmax6(&a8,&a9,&a10,&a11,&a12,&a13);
3505   minmax6(&b8,&b9,&b10,&b11,&b12,&b13);
3506
3507   //chargement valeur suivante (23)
3508   a13 = tex2D(tex_img_inc, j-2, i  );
3509   b13 = tex2D(tex_img_inc, j+3, i  );
3510   //minmax aux extremites
3511   minmax5(&a9,&a10,&a11,&a12,&a13);
3512   minmax5(&b9,&b10,&b11,&b12,&b13);
3513
3514   //chargement valeur suivante (24)
3515   a13 = tex2D(tex_img_inc, j-2, i+1);
3516   b13 = tex2D(tex_img_inc, j+3, i+1);
3517   //minmax aux extremites
3518   minmax4(&a10,&a11,&a12,&a13);
3519   minmax4(&b10,&b11,&b12,&b13);
3520
3521   //chargement valeur suivante (25)
3522   a13 = tex2D(tex_img_inc, j-2, i+2);
3523   b13 = tex2D(tex_img_inc, j+3, i+2);
3524   //minmax aux extremites
3525   minmax3(&a11,&a12,&a13);
3526   minmax3(&b11,&b12,&b13);
3527
3528   //median au milieu !
3529   output[ __mul24(i, j_dim) +j ] = a12 ;
3530   output[ __mul24(i, j_dim) +j+1 ] = b12 ;
3531
3532 }
3533
3534
3535 __global__ void kernel_median7( unsigned short*output, int i_dim, int j_dim)
3536 {
3537
3538   // coordonnees absolues du point
3539   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
3540   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
3541
3542   /**************************************************************************
3543    *             tri(s)
3544    **************************************************************************/
3545   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 ;
3546
3547   /****************************************
3548    * les 26 premieres valeurs
3549    ****************************************/
3550   //premiere ligne
3551   a0  = tex2D(tex_img_ins, j-2, i-3) ;
3552   a1  = tex2D(tex_img_ins, j-1, i-3) ;
3553   a2  = tex2D(tex_img_ins, j  , i-3) ;
3554   a3  = tex2D(tex_img_ins, j+1, i-3) ;
3555   a4  = tex2D(tex_img_ins, j+2, i-3) ;
3556   a5  = tex2D(tex_img_ins, j+3, i-3) ;
3557   //deuxieme ligne
3558   a6  = tex2D(tex_img_ins, j-2, i-2) ;
3559   a7  = tex2D(tex_img_ins, j-1, i-2) ;
3560   a8  = tex2D(tex_img_ins, j  , i-2) ;
3561   a9  = tex2D(tex_img_ins, j+1, i-2) ;
3562   a10 = tex2D(tex_img_ins, j+2, i-2) ;
3563   a11 = tex2D(tex_img_ins, j+3, i-2) ;
3564   //troisieme ligne
3565   a12 = tex2D(tex_img_ins, j-2, i-1) ;
3566   a13 = tex2D(tex_img_ins, j-1, i-1) ;
3567   a14 = tex2D(tex_img_ins, j  , i-1) ;
3568   a15 = tex2D(tex_img_ins, j+1, i-1) ;
3569   a16 = tex2D(tex_img_ins, j+2, i-1) ;
3570   a17 = tex2D(tex_img_ins, j+3, i-1) ;
3571   //quatrieme ligne
3572   a18 = tex2D(tex_img_ins, j-2, i  ) ;
3573   a19 = tex2D(tex_img_ins, j-1, i  ) ;
3574   a20 = tex2D(tex_img_ins, j  , i  ) ;
3575   a21 = tex2D(tex_img_ins, j+1, i  ) ;
3576   a22 = tex2D(tex_img_ins, j+2, i  ) ;
3577   a23 = tex2D(tex_img_ins, j+3, i  ) ;
3578   //cinquieme ligne
3579   a24 = tex2D(tex_img_ins, j-2, i+1) ;
3580   a25 = tex2D(tex_img_ins, j-1, i+1) ;
3581
3582   //min max aux extremites
3583   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);
3584
3585   //chargement valeurs suivante (26)
3586   a25 = tex2D(tex_img_ins, j , i+1);
3587   //minmax aux extremites
3588   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);
3589
3590   //chargement valeur suivante (27)
3591   a25 = tex2D(tex_img_ins, j+1, i+1);
3592   //minmax aux extremites
3593   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);
3594
3595   //chargement valeur suivante (28)
3596   a25 = tex2D(tex_img_ins, j+2, i+1);
3597   //minmax aux extremites
3598   minmax23(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3599
3600   //chargement valeur suivante (29)
3601   a25 = tex2D(tex_img_ins, j+3, i+1);
3602   //minmax aux extremites
3603   minmax22(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3604
3605   //chargement valeur suivante (30)
3606   a25 = tex2D(tex_img_ins, j-2, i+2);
3607   //minmax aux extremites
3608   minmax21(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3609
3610   //chargement valeur suivante (31)
3611   a25 = tex2D(tex_img_ins, j-1, i+2);
3612   //minmax aux extmites
3613   minmax20(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3614
3615   //chargement valeur suivante (32)
3616   a25 = tex2D(tex_img_ins, j  , i+2);
3617   //minmax aux extmites
3618   minmax19(&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3619
3620   //chargement valeur suivante (33)
3621   a25 = tex2D(tex_img_ins, j+1, i+2);
3622   //minmax aux extmites
3623   minmax18(&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3624
3625   //chargement valeur suivante (34)
3626   a25 = tex2D(tex_img_ins, j+2, i+2);
3627   //minmax aux extmites
3628   minmax17(&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3629
3630   //chargement valeur suivante (35)
3631   a25 = tex2D(tex_img_ins, j+3, i+2);
3632   //minmax aux extmites
3633   minmax16(&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3634
3635   //chargement valeur suivante (36)
3636   a25 = tex2D(tex_img_ins, j-2, i+3);
3637   //minmax aux extmites
3638   minmax15(&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3639
3640   //chargement valeur suivante (37)
3641   a25 = tex2D(tex_img_ins, j-1, i+3);
3642   //minmax aux extmites
3643
3644   minmax14(&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3645   //chargement valeur suivante (38)
3646   a25 = tex2D(tex_img_ins, j  , i+3);
3647   //minmax aux extmites
3648   minmax13(&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3649
3650   //chargement valeur suivante (39)
3651   a25 = tex2D(tex_img_ins, j+1, i+3);
3652   //minmax aux extmites
3653   minmax12(&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3654
3655   //chargement valeur suivante (40)
3656   a25 = tex2D(tex_img_ins, j+2, i+3);
3657   //minmax aux extmites
3658   minmax11(&a15, &a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3659
3660   //chargement valeur suivante (41)
3661   a25 = tex2D(tex_img_ins, j+3, i+3);
3662   //minmax aux extmites
3663   minmax10(&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3664
3665   //chargement valeurs suivantes
3666   a25 = tex2D(tex_img_ins, j-3, i+3);
3667   //minmax aux extremites
3668   minmax9(&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3669
3670   //chargement valeurs suivantes
3671   a25 = tex2D(tex_img_ins, j-3, i+2);
3672   //minmax aux extremites
3673   minmax8(&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3674
3675   //chargement valeurs suivantes
3676   a25 = tex2D(tex_img_ins, j-3, i+1);
3677   //minmax aux extremites
3678   minmax7(&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3679
3680   //chargement valeurs suivantes
3681   a25 = tex2D(tex_img_ins, j-3, i);
3682   //minmax aux extremites
3683   minmax6(&a20,&a21,&a22,&a23,&a24,&a25);
3684
3685   //chargement valeurs suivantes
3686   a25 = tex2D(tex_img_ins, j-3, i-1);
3687   //minmax aux extremites
3688   minmax5(&a21,&a22,&a23,&a24,&a25);
3689
3690   //chargement valeurs suivantes
3691   a25 = tex2D(tex_img_ins, j-3, i-2);
3692   //minmax aux extremites
3693   minmax4(&a22,&a23,&a24,&a25);
3694
3695   //chargement valeurs suivantes
3696   a25 = tex2D(tex_img_ins, j-3, i-3);
3697   //minmax aux extremites
3698   minmax3(&a23,&a24,&a25);
3699
3700   //medians au milieu !
3701   output[ __mul24(i, j_dim) +j ] = a24 ;
3702
3703 }
3704
3705
3706
3707 __global__ void kernel_median7_2pix( unsigned short*output, int i_dim, int j_dim)
3708 {
3709
3710   // coordonnees absolues du point
3711   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
3712   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
3713
3714   /**************************************************************************
3715    *             tri(s)
3716    **************************************************************************/
3717   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 ;
3718   int b17,b18,b19,b20,b21,b22,b23,b24,b25;
3719
3720   /****************************************
3721    * les 26 premieres valeurs
3722    ****************************************/
3723   //premiere ligne
3724   a0  = tex2D(tex_img_ins, j-2, i-3) ;
3725   a1  = tex2D(tex_img_ins, j-1, i-3) ;
3726   a2  = tex2D(tex_img_ins, j  , i-3) ;
3727   a3  = tex2D(tex_img_ins, j+1, i-3) ;
3728   a4  = tex2D(tex_img_ins, j+2, i-3) ;
3729   a5  = tex2D(tex_img_ins, j+3, i-3) ;
3730   //deuxieme ligne
3731   a6  = tex2D(tex_img_ins, j-2, i-2) ;
3732   a7  = tex2D(tex_img_ins, j-1, i-2) ;
3733   a8  = tex2D(tex_img_ins, j  , i-2) ;
3734   a9  = tex2D(tex_img_ins, j+1, i-2) ;
3735   a10 = tex2D(tex_img_ins, j+2, i-2) ;
3736   a11 = tex2D(tex_img_ins, j+3, i-2) ;
3737   //troisieme ligne
3738   a12 = tex2D(tex_img_ins, j-2, i-1) ;
3739   a13 = tex2D(tex_img_ins, j-1, i-1) ;
3740   a14 = tex2D(tex_img_ins, j  , i-1) ;
3741   a15 = tex2D(tex_img_ins, j+1, i-1) ;
3742   a16 = tex2D(tex_img_ins, j+2, i-1) ;
3743   a17 = tex2D(tex_img_ins, j+3, i-1) ;
3744   //quatrieme ligne
3745   a18 = tex2D(tex_img_ins, j-2, i  ) ;
3746   a19 = tex2D(tex_img_ins, j-1, i  ) ;
3747   a20 = tex2D(tex_img_ins, j  , i  ) ;
3748   a21 = tex2D(tex_img_ins, j+1, i  ) ;
3749   a22 = tex2D(tex_img_ins, j+2, i  ) ;
3750   a23 = tex2D(tex_img_ins, j+3, i  ) ;
3751   //cinquieme ligne
3752   a24 = tex2D(tex_img_ins, j-2, i+1) ;
3753   a25 = tex2D(tex_img_ins, j-1, i+1) ;
3754
3755   //min max aux extremites
3756   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);
3757
3758   //chargement valeurs suivante (26)
3759   a25 = tex2D(tex_img_ins, j , i+1);
3760   //minmax aux extremites
3761   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);
3762
3763   //chargement valeur suivante (27)
3764   a25 = tex2D(tex_img_ins, j+1, i+1);
3765   //minmax aux extremites
3766   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);
3767
3768   //chargement valeur suivante (28)
3769   a25 = tex2D(tex_img_ins, j+2, i+1);
3770   //minmax aux extremites
3771   minmax23(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3772
3773   //chargement valeur suivante (29)
3774   a25 = tex2D(tex_img_ins, j+3, i+1);
3775   //minmax aux extremites
3776   minmax22(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3777
3778   //chargement valeur suivante (30)
3779   a25 = tex2D(tex_img_ins, j-2, i+2);
3780   //minmax aux extremites
3781   minmax21(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3782
3783   //chargement valeur suivante (31)
3784   a25 = tex2D(tex_img_ins, j-1, i+2);
3785   //minmax aux extmites
3786   minmax20(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3787
3788   //chargement valeur suivante (32)
3789   a25 = tex2D(tex_img_ins, j  , i+2);
3790   //minmax aux extmites
3791   minmax19(&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3792
3793   //chargement valeur suivante (33)
3794   a25 = tex2D(tex_img_ins, j+1, i+2);
3795   //minmax aux extmites
3796   minmax18(&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3797
3798   //chargement valeur suivante (34)
3799   a25 = tex2D(tex_img_ins, j+2, i+2);
3800   //minmax aux extmites
3801   minmax17(&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3802
3803   //chargement valeur suivante (35)
3804   a25 = tex2D(tex_img_ins, j+3, i+2);
3805   //minmax aux extmites
3806   minmax16(&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3807
3808   //chargement valeur suivante (36)
3809   a25 = tex2D(tex_img_ins, j-2, i+3);
3810   //minmax aux extmites
3811   minmax15(&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3812
3813   //chargement valeur suivante (37)
3814   a25 = tex2D(tex_img_ins, j-1, i+3);
3815   //minmax aux extmites
3816
3817   minmax14(&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3818   //chargement valeur suivante (38)
3819   a25 = tex2D(tex_img_ins, j  , i+3);
3820   //minmax aux extmites
3821   minmax13(&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3822
3823   //chargement valeur suivante (39)
3824   a25 = tex2D(tex_img_ins, j+1, i+3);
3825   //minmax aux extmites
3826   minmax12(&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3827
3828   //chargement valeur suivante (40)
3829   a25 = tex2D(tex_img_ins, j+2, i+3);
3830   //minmax aux extmites
3831   minmax11(&a15, &a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3832
3833   //chargement valeur suivante (41)
3834   a25 = tex2D(tex_img_ins, j+3, i+3);
3835   //minmax aux extmites
3836   minmax10(&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3837
3838   // fin des pixels voisins communs deux pixels centraux
3839   b17=a17; b18=a18; b19=a19; b20=a20; b21=a21; b22=a22; b23=a23; b24=a24; b25=a25;
3840
3841   //chargement valeurs suivantes
3842   a25 = tex2D(tex_img_ins, j-3, i+3);
3843   b25 = tex2D(tex_img_ins, j+4, i+3);
3844   //minmax aux extremites
3845   minmax9(&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3846   minmax9(&b17,&b18,&b19,&b20,&b21,&b22,&b23,&b24,&b25);
3847
3848   //chargement valeurs suivantes
3849   a25 = tex2D(tex_img_ins, j-3, i+2);
3850   b25 = tex2D(tex_img_ins, j+4, i+2);
3851   //minmax aux extremites
3852   minmax8(&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3853   minmax8(&b18,&b19,&b20,&b21,&b22,&b23,&b24,&b25);
3854
3855   //chargement valeurs suivantes
3856   a25 = tex2D(tex_img_ins, j-3, i+1);
3857   b25 = tex2D(tex_img_ins, j+4, i+1);
3858   //minmax aux extremites
3859   minmax7(&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3860   minmax7(&b19,&b20,&b21,&b22,&b23,&b24,&b25);
3861
3862   //chargement valeurs suivantes
3863   a25 = tex2D(tex_img_ins, j-3, i);
3864   b25 = tex2D(tex_img_ins, j+4, i);
3865   //minmax aux extremites
3866   minmax6(&a20,&a21,&a22,&a23,&a24,&a25);
3867   minmax6(&b20,&b21,&b22,&b23,&b24,&b25);
3868
3869   //chargement valeurs suivantes
3870   a25 = tex2D(tex_img_ins, j-3, i-1);
3871   b25 = tex2D(tex_img_ins, j+4, i-1);
3872   //minmax aux extremites
3873   minmax5(&a21,&a22,&a23,&a24,&a25);
3874   minmax5(&b21,&b22,&b23,&b24,&b25);
3875
3876   //chargement valeurs suivantes
3877   a25 = tex2D(tex_img_ins, j-3, i-2);
3878   b25 = tex2D(tex_img_ins, j+4, i-2);
3879   //minmax aux extremites
3880   minmax4(&a22,&a23,&a24,&a25);
3881   minmax4(&b22,&b23,&b24,&b25);
3882
3883   //chargement valeurs suivantes
3884   a25 = tex2D(tex_img_ins, j-3, i-3);
3885   b25 = tex2D(tex_img_ins, j+4, i-3);
3886   //minmax aux extremites
3887   minmax3(&a23,&a24,&a25);
3888   minmax3(&b23,&b24,&b25);
3889
3890   //medians au milieu !
3891   output[ __mul24(i, j_dim) +j ] = a24 ;
3892   output[ __mul24(i, j_dim) +j+1 ] = b24 ;
3893
3894 }
3895
3896
3897 __global__ void kernel_median7_2pix( unsigned char  *output, int i_dim, int j_dim)
3898 {
3899
3900   // coordonnees absolues du point
3901   int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ;
3902   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
3903
3904   /**************************************************************************
3905    *             tri(s)
3906    **************************************************************************/
3907   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 ;
3908   int b17,b18,b19,b20,b21,b22,b23,b24,b25;
3909
3910   /****************************************
3911    * les 26 premieres valeurs
3912    ****************************************/
3913   //premiere ligne
3914   a0  = tex2D(tex_img_inc, j-2, i-3) ;
3915   a1  = tex2D(tex_img_inc, j-1, i-3) ;
3916   a2  = tex2D(tex_img_inc, j  , i-3) ;
3917   a3  = tex2D(tex_img_inc, j+1, i-3) ;
3918   a4  = tex2D(tex_img_inc, j+2, i-3) ;
3919   a5  = tex2D(tex_img_inc, j+3, i-3) ;
3920   //deuxieme ligne
3921   a6  = tex2D(tex_img_inc, j-2, i-2) ;
3922   a7  = tex2D(tex_img_inc, j-1, i-2) ;
3923   a8  = tex2D(tex_img_inc, j  , i-2) ;
3924   a9  = tex2D(tex_img_inc, j+1, i-2) ;
3925   a10 = tex2D(tex_img_inc, j+2, i-2) ;
3926   a11 = tex2D(tex_img_inc, j+3, i-2) ;
3927   //troisieme ligne
3928   a12 = tex2D(tex_img_inc, j-2, i-1) ;
3929   a13 = tex2D(tex_img_inc, j-1, i-1) ;
3930   a14 = tex2D(tex_img_inc, j  , i-1) ;
3931   a15 = tex2D(tex_img_inc, j+1, i-1) ;
3932   a16 = tex2D(tex_img_inc, j+2, i-1) ;
3933   a17 = tex2D(tex_img_inc, j+3, i-1) ;
3934   //quatrieme ligne
3935   a18 = tex2D(tex_img_inc, j-2, i  ) ;
3936   a19 = tex2D(tex_img_inc, j-1, i  ) ;
3937   a20 = tex2D(tex_img_inc, j  , i  ) ;
3938   a21 = tex2D(tex_img_inc, j+1, i  ) ;
3939   a22 = tex2D(tex_img_inc, j+2, i  ) ;
3940   a23 = tex2D(tex_img_inc, j+3, i  ) ;
3941   //cinquieme ligne
3942   a24 = tex2D(tex_img_inc, j-2, i+1) ;
3943   a25 = tex2D(tex_img_inc, j-1, i+1) ;
3944
3945   //min max aux extremites
3946   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);
3947
3948   //chargement valeurs suivante (26)
3949   a25 = tex2D(tex_img_inc, j , i+1);
3950   //minmax aux extremites
3951   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);
3952
3953   //chargement valeur suivante (27)
3954   a25 = tex2D(tex_img_inc, j+1, i+1);
3955   //minmax aux extremites
3956   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);
3957
3958   //chargement valeur suivante (28)
3959   a25 = tex2D(tex_img_inc, j+2, i+1);
3960   //minmax aux extremites
3961   minmax23(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3962
3963   //chargement valeur suivante (29)
3964   a25 = tex2D(tex_img_inc, j+3, i+1);
3965   //minmax aux extremites
3966   minmax22(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3967
3968   //chargement valeur suivante (30)
3969   a25 = tex2D(tex_img_inc, j-2, i+2);
3970   //minmax aux extremites
3971   minmax21(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3972
3973   //chargement valeur suivante (31)
3974   a25 = tex2D(tex_img_inc, j-1, i+2);
3975   //minmax aux extmites
3976   minmax20(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3977
3978   //chargement valeur suivante (32)
3979   a25 = tex2D(tex_img_inc, j  , i+2);
3980   //minmax aux extmites
3981   minmax19(&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3982
3983   //chargement valeur suivante (33)
3984   a25 = tex2D(tex_img_inc, j+1, i+2);
3985   //minmax aux extmites
3986   minmax18(&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3987
3988   //chargement valeur suivante (34)
3989   a25 = tex2D(tex_img_inc, j+2, i+2);
3990   //minmax aux extmites
3991   minmax17(&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3992
3993   //chargement valeur suivante (35)
3994   a25 = tex2D(tex_img_inc, j+3, i+2);
3995   //minmax aux extmites
3996   minmax16(&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
3997
3998   //chargement valeur suivante (36)
3999   a25 = tex2D(tex_img_inc, j-2, i+3);
4000   //minmax aux extmites
4001   minmax15(&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
4002
4003   //chargement valeur suivante (37)
4004   a25 = tex2D(tex_img_inc, j-1, i+3);
4005   //minmax aux extmites
4006
4007   minmax14(&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
4008   //chargement valeur suivante (38)
4009   a25 = tex2D(tex_img_inc, j  , i+3);
4010   //minmax aux extmites
4011   minmax13(&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
4012
4013   //chargement valeur suivante (39)
4014   a25 = tex2D(tex_img_inc, j+1, i+3);
4015   //minmax aux extmites
4016   minmax12(&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
4017
4018   //chargement valeur suivante (40)
4019   a25 = tex2D(tex_img_inc, j+2, i+3);
4020   //minmax aux extmites
4021   minmax11(&a15, &a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
4022
4023   //chargement valeur suivante (41)
4024   a25 = tex2D(tex_img_inc, j+3, i+3);
4025   //minmax aux extmites
4026   minmax10(&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
4027
4028   // fin des pixels voisins communs deux pixels centraux
4029   b17=a17; b18=a18; b19=a19; b20=a20; b21=a21; b22=a22; b23=a23; b24=a24; b25=a25;
4030
4031   //chargement valeurs suivantes
4032   a25 = tex2D(tex_img_inc, j-3, i+3);
4033   b25 = tex2D(tex_img_inc, j+4, i+3);
4034   //minmax aux extremites
4035   minmax9(&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
4036   minmax9(&b17,&b18,&b19,&b20,&b21,&b22,&b23,&b24,&b25);
4037
4038   //chargement valeurs suivantes
4039   a25 = tex2D(tex_img_inc, j-3, i+2);
4040   b25 = tex2D(tex_img_inc, j+4, i+2);
4041   //minmax aux extremites
4042   minmax8(&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
4043   minmax8(&b18,&b19,&b20,&b21,&b22,&b23,&b24,&b25);
4044
4045   //chargement valeurs suivantes
4046   a25 = tex2D(tex_img_inc, j-3, i+1);
4047   b25 = tex2D(tex_img_inc, j+4, i+1);
4048   //minmax aux extremites
4049   minmax7(&a19,&a20,&a21,&a22,&a23,&a24,&a25);
4050   minmax7(&b19,&b20,&b21,&b22,&b23,&b24,&b25);
4051
4052   //chargement valeurs suivantes
4053   a25 = tex2D(tex_img_inc, j-3, i);
4054   b25 = tex2D(tex_img_inc, j+4, i);
4055   //minmax aux extremites
4056   minmax6(&a20,&a21,&a22,&a23,&a24,&a25);
4057   minmax6(&b20,&b21,&b22,&b23,&b24,&b25);
4058
4059   //chargement valeurs suivantes
4060   a25 = tex2D(tex_img_inc, j-3, i-1);
4061   b25 = tex2D(tex_img_inc, j+4, i-1);
4062   //minmax aux extremites
4063   minmax5(&a21,&a22,&a23,&a24,&a25);
4064   minmax5(&b21,&b22,&b23,&b24,&b25);
4065
4066   //chargement valeurs suivantes
4067   a25 = tex2D(tex_img_inc, j-3, i-2);
4068   b25 = tex2D(tex_img_inc, j+4, i-2);
4069   //minmax aux extremites
4070   minmax4(&a22,&a23,&a24,&a25);
4071   minmax4(&b22,&b23,&b24,&b25);
4072
4073   //chargement valeurs suivantes
4074   a25 = tex2D(tex_img_inc, j-3, i-3);
4075   b25 = tex2D(tex_img_inc, j+4, i-3);
4076   //minmax aux extremites
4077   minmax3(&a23,&a24,&a25);
4078   minmax3(&b23,&b24,&b25);
4079
4080   //medians au milieu !
4081   output[ __mul24(i, j_dim) +j ] = a24 ;
4082   output[ __mul24(i, j_dim) +j+1 ] = b24 ;
4083
4084 }
4085
4086
4087
4088 /*****************************************************************************
4089  *      median generic shared mem + forgetfull
4090  *****************************************************************************/
4091 __global__ void kernel_medianForgetRSH( unsigned short * output, int i_dim, int j_dim, int r)
4092 {
4093   // coordonnees absolues du point
4094   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
4095   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
4096   unsigned short ic, jc ;
4097   unsigned short Nreg = ((2*r+1)*(2*r+1))/2 + 2 ;
4098
4099   int bdimX = blockDim.x; //<<1 ; // pour le cas \E0 2 pixels par thread
4100   int tidX = threadIdx.x; //<<1 ;
4101
4102   // chargement en smem
4103   int idrow = threadIdx.y*(blockDim.x+2*r) ;
4104
4105   extern __shared__ int medRoi[];
4106
4107   // bloc 0 (en haut \E0 gauche)
4108   medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-r, i-r) ;
4109   // bloc 1 (en haut \E0 droite)...
4110   if ( threadIdx.x < 2*r ) //...ou plutot ce qu'il en manque
4111         medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-r, i-r ) ;
4112   // bloc 2 ( en bas \E0 gauche)
4113   if ( threadIdx.y < 2*r )
4114         {
4115           idrow = (threadIdx.y+blockDim.y)*(blockDim.x+2*r) ;
4116           medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-r, i+blockDim.y-r ) ;
4117           //bloc 4 ( en bas \E0 doite )...
4118           if ( threadIdx.x < 2*r ) //...ou ce qu'il en manque
4119                 medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-r, i+blockDim.y-r ) ;
4120         }
4121   __syncthreads() ;
4122
4123   // remplissage du vecteur de tri minmax
4124   unsigned short vect[8066] ;
4125   int Freg=Nreg ;
4126   for (ic=0; ic<2*r+1; ic++)
4127         {
4128           for (jc=0; jc<2*r+1; jc++)
4129                 {
4130                   if ( ic*(2*r+1)+jc < Nreg )
4131                         {
4132                           vect[ ic*(2*r+1)+jc ] = medRoi[ (threadIdx.y+ic)*(bdimX+2*r)+ (tidX+jc) ] ;
4133                         } else
4134                         {
4135                           minmaxN(vect, Freg--) ;
4136                           vect[ Nreg-1 ] = medRoi[ (threadIdx.y+ic)*(bdimX+2*r)+ (tidX+jc) ] ;
4137                         }
4138                 }
4139         }
4140   minmax3(&vect[Nreg-3], &vect[Nreg-2], &vect[Nreg-1])
4141
4142   //medRoi[ (threadIdx.y+ic)*(bdimX+L-1)+ (tidX+jc) ]
4143
4144         output[ __mul24(i, j_dim) +j ] = vect[ Nreg-2 ];
4145 }
4146
4147
4148 /*****************************************************************************
4149  *      median generic shared mem + forgetfull
4150  *****************************************************************************/
4151 __global__ void kernel_medianForgetR( unsigned short * output, int i_dim, int j_dim, int r)
4152 {
4153   // coordonnees absolues du point
4154   int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ;
4155   int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
4156   unsigned short ic, jc ;
4157   unsigned short Nreg = ((2*r+1)*(2*r+1))/2 + 2 ;
4158
4159   // remplissage du vecteur de tri minmax
4160   unsigned short vect[8066] ;
4161   int Freg=Nreg ;
4162   for (ic=0; ic<2*r+1; ic++)
4163         {
4164           for (jc=0; jc<2*r+1; jc++)
4165                 {
4166                   if ( ic*(2*r+1)+jc < Nreg )
4167                         {
4168                           vect[ ic*(2*r+1)+jc ] = tex2D(tex_img_ins, j-r+jc, i-r+ic) ;
4169                         } else
4170                         {
4171                           minmaxN(vect, Freg--) ;
4172                           vect[ Nreg-1 ] = tex2D(tex_img_ins, j-r+jc, i-r+ic) ;
4173                         }
4174                 }
4175         }
4176   minmax3(&vect[Nreg-3], &vect[Nreg-2], &vect[Nreg-1])
4177
4178   //medRoi[ (threadIdx.y+ic)*(bdimX+L-1)+ (tidX+jc) ]
4179
4180         output[ __mul24(i, j_dim) +j ] = vect[ Nreg-2 ];
4181 }
4182
4183
4184
4185 /**********************************************************************************
4186  *    MEDIAN PSEUDO-SEPARABLE POUR GRANDS NOYAUX
4187  **********************************************************************************/
4188 __global__ void kernel_medianSepR( unsigned short *output, int i_dim, int j_dim, int r)
4189 {
4190
4191   int idc, val, min, max, inf, egal, sup, mxinf, minsup, estim ;
4192
4193   //coordonnées ds le bloc
4194   int ib = threadIdx.y ;
4195   int jb = threadIdx.x ;
4196   //int idx_h = __mul24(ib+r,blockDim.x) + jb ;   // index pixel deans shmem (bloc+halo)
4197   //int offset = __mul24(blockDim.x,r) ;
4198
4199   // coordonnees absolues du point
4200   int j = __mul24(blockIdx.x,blockDim.x) + jb ;
4201   int i = __mul24(blockIdx.y,blockDim.y) + ib ;
4202
4203   extern __shared__ int buff[] ;
4204   /***********************************************************************************
4205    *              CHARGEMENT DATA EN SHARED MEM
4206    ***********************************************************************************/
4207   for (idc = 0 ; idc < (2*(blockDim.y+r)-1)/blockDim.y; idc++)
4208         {
4209           if (idc*blockDim.y +ib < i_dim)
4210                 buff[ (idc*blockDim.y +ib)*blockDim.x + jb ]    = tex2D(tex_img_ins, j, i-r+ idc*blockDim.y) ;
4211         }
4212
4213   __syncthreads() ;
4214   /**********************************************************************************************
4215    *               TRI VERTICAL par algo TORBEN MOGENSEN
4216    *          (a little bit slow but saves memory => faster !)
4217    **********************************************************************************************/
4218   min = max = buff[ ib*blockDim.x +jb] ;
4219
4220   for (idc= 0 ; idc< 2*r+1 ; idc++ )
4221         {
4222           val = buff[ __mul24(ib+idc, blockDim.x) +jb ] ;
4223           if ( val  < min ) min = val ;
4224           if ( val  > max ) max = val ;
4225         }
4226
4227   while (1)
4228         {
4229           estim = (min+max)/2 ;
4230           inf = sup = egal = 0  ;
4231           mxinf = min ;
4232           minsup= max ;
4233           for (idc =0; idc< 2*r+1 ; idc++)
4234                 {
4235                   val = buff[ __mul24(ib+idc, blockDim.x) +jb ] ;
4236                   if( val < estim )
4237                         {
4238                           inf++;
4239                           if( val > mxinf) mxinf = val ;
4240                         } else if (val > estim)
4241                         {
4242                           sup++;
4243                           if( val < minsup) minsup = val ;
4244                         } else egal++ ;
4245                 }
4246           if ( (inf <= (r+1))&&(sup <=(r+1)) ) break ;
4247           else if (inf>sup) max = mxinf ;
4248           else min = minsup ;
4249           }
4250
4251   if ( inf >= r+1 ) val = mxinf ;
4252   else if (inf+egal >= r+1) val = estim ;
4253   else val = minsup ;
4254
4255   output[ __mul24(j, i_dim)  +i  ] =  val ;
4256 }
4257
4258
4259 /**
4260  *
4261  * correction de staircase
4262  */
4263 __global__ void kernel_staircase_reduc3(unsigned int * img_out, unsigned int L, unsigned int H)
4264 {
4265   // coordonnees du point dans le bloc
4266   //unsigned int iib = threadIdx.x ;
4267   //unsigned int jib = threadIdx.y ;
4268   // coordonnees du point dans l'image
4269   unsigned int y = blockIdx.y*blockDim.y + threadIdx.y;
4270   unsigned int x = blockIdx.x*blockDim.x + threadIdx.x;
4271
4272   int a, b, c, d, e, f, g, h, i ;            // gl des voisins
4273   float wa, wb, wc, wd, we, wf, wg, wh, wi ; // poids
4274   float S1, S2, S11, S22, S12, S0, Sx, Sx1, Sx2 ;
4275   float c1,c2 ;
4276
4277   // chargement des valeurs GL des pixels du voisinage
4278   a = tex2D(tex_img_estim, x-1, y-1) ;
4279   b = tex2D(tex_img_estim, x  , y-1) ;
4280   c = tex2D(tex_img_estim, x+1, y-1) ;
4281   d = tex2D(tex_img_estim, x-1, y  ) ;
4282   e = tex2D(tex_img_estim, x  , y  ) ;
4283   f = tex2D(tex_img_estim, x+1, y  ) ;
4284   g = tex2D(tex_img_estim, x-1, y+1) ;
4285   h = tex2D(tex_img_estim, x  , y+1) ;
4286   i = tex2D(tex_img_estim, x+1, y+1) ;
4287
4288   wa = tex1D(tex_noyau, abs(a-e)) ;
4289   wb = tex1D(tex_noyau, abs(b-e)) ;
4290   wc = tex1D(tex_noyau, abs(c-e)) ;
4291   wd = tex1D(tex_noyau, abs(d-e)) ;
4292   we = 1 ;
4293   wf = tex1D(tex_noyau, abs(f-e)) ;
4294   wg = tex1D(tex_noyau, abs(g-e)) ;
4295   wh = tex1D(tex_noyau, abs(h-e)) ;
4296   wi = tex1D(tex_noyau, abs(i-e)) ;
4297
4298
4299   //construction des elements du systeme lineaire
4300   S0  = wa+wb+wc+wd+we+wf+wg+wh+wi ;
4301   S1  = wc+wf+wi-wa-wd-wg ;
4302   S2  = wg+wh+wi-wa-wb-wc ;
4303   Sx  = wa*a + wb*b + wc*c + wd*d + we*e + wf*f + wg*g + wh*h + wi*i ;
4304   Sx1 = wc*c + wf*f + wi*i - wa*a - wd*d - wg*g ;
4305   Sx2 = wg*g + wh*h + wi*i - wa*a - wb*b - wc*c ;
4306   S11 = wc+wf+wi+wa+wd+wg ;
4307   S22 = wg+wh+wi+wa+wb+wc ;
4308   S12 = wa + wi -wc -wg ;
4309
4310   c1 = S22*(S11*Sx-S1*Sx1) + S12*(S1*Sx2-S12*Sx) + S2*(S12*Sx1-S11*Sx2) ;
4311   c2 = S22*(S11*S0-S1*S1) + S12*(S2*S1-S12*S0) + S2*(S1*S12-S2*S11) ;
4312   img_out[y*L + x] = c1/c2 ;
4313 }
4314
4315