levelines_kernels.cu

   1 /************************************************************************
   2  * chemins des lignes de niveaux pour la version à chemins constants
   3  * Ne sont conservés que pour comparaison GPU/CPU -> à faire disparaître
   4  * longueur = 4 pixels
   5  * une ligne = un chemin
   6  ************************************************************************/
   7 __constant__ int pathDi[PSIZE_I][PSIZE_J-1] =
   8   {
   9     // Q1
  10         {  -1, -1, -1},       // 90
  11         {  -1, -1, -1},       // 75
  12         {  -1, -1, -1},       // 60
  13         {  -1, -1, -1},       // 45
  14         {  -1,  0, -1},       // 30
  15         {   0, -1,  0},       // 15
  16         // Q4
  17         {   0,  0,  0},       // 0
  18         {   0,  1,  0},       // 345
  19         {   1,  0,  1},       // 330
  20         {   1,  1,  1},       // 315
  21         {   1,  1,  1},       // 300
  22         {   1,  1,  1},       // 285
  23         // Q3
  24         {   1,  1,  1},       // 270
  25         {   1,  1,  1},       // 255
  26         {   1,  1,  1},       // 240
  27         {   1,  1,  1},       // 225
  28         {   1,  0,  1},       // 210
  29         {   0,  1,  0},       // 195
  30         // Q2
  31         {   0,  0,  0},       // 180
  32         {   0, -1,  0},       // 165
  33         {  -1,  0, -1},       // 150
  34         {  -1, -1, -1},       // 135
  35         {  -1, -1, -1},       // 120
  36         {  -1, -1, -1}        // 105
  37   } ;     //
  38
  39 __constant__ int pathDj[PSIZE_I][PSIZE_J-1] =
  40   {
  41         // Q1
  42         {  0,  0,  0},       // 90
  43         {  0,  1,  0},       // 75
  44         {  1,  0,  1},       // 60
  45         {  1,  1,  1},       // 45
  46         {  1,  1,  1},       // 30
  47         {  1,  1,  1},       // 15
  48         // Q4
  49         {  1,  1,  1},       // 0
  50         {  1,  1,  1},       // 345
  51         {  1,  1,  1},       // 330
  52         {  1,  1,  1},       // 315
  53         {  1,  0,  1},       // 300
  54         {  0,  1,  0},       // 285
  55         // Q3
  56         {  0,  0,  0},       // 270
  57         {  0, -1,  0},       // 255
  58         { -1,  0, -1},       // 240
  59         { -1, -1, -1},       // 225
  60         { -1, -1, -1},       // 210
  61         { -1, -1, -1},       // 195
  62         // Q2
  63         { -1, -1, -1},       // 180
  64         { -1, -1, -1},       // 165
  65         { -1, -1, -1},       // 150
  66         { -1, -1, -1},       // 135
  67         { -1,  0, -1},       // 120
  68         {  0, -1,  0}        // 105
  69   } ;
  70
  71
  72
  73 // valeurs des tangentes des angles de base pour la génération initiale des chemins
  74 // pour la version à chemins de longueur paramétrable
  75 __constant__ float tangente[] = {0.000, 0.268, 0.577, 1.000} ;
  76
  77 // declarations des textures
  78 texture<int, 2, cudaReadModeElementType> tex_img_in ;
  79 texture<int, 2, cudaReadModeElementType> tex_img_estim ;
  80 texture<int, 2, cudaReadModeElementType> tex_img_lniv ;
  81 texture<int2, 2, cudaReadModeElementType> tex_paths ;
  82
  83
  84
  85 /**
  86  *
  87  * \brief calcule les chemins
  88  * \author NB - PhyTI, modifié by zulu pour adaptater aux chemins paramétrables
  89  *
  90  * \param[in] r         longueur des chemins
  91  *
  92  * \param[out] d_paths  matrice des déplacements relatifs (chemins)
  93  *
  94  * Cette fonction utilise le tableau constant des tangentes des angles
  95  * considérés pour le calcul de chemins (float tangente[]).
  96  *
  97  */
  98 __global__ void kernel_calcul_paths( int2 * d_paths, unsigned int r){
  99
 100   unsigned int idpath = 0 ;
 101   int ic, jc, iprec, jprec ;
 102   float offset = 0.5 ;
 103   unsigned int basepath = 0 ;
 104
 105   // Q1 inf
 106   for (int a=0 ; a< 4 ; a++){        // les 4 angles 0,15,30 et 45
 107         for (int p=0 ; p< r ; p++){      // les r points
 108           ic = r-1 - floor(tangente[a]*p + offset) ;
 109           if ( p > 0 ){
 110                 d_paths[idpath*(r-1)+p-1].x = ic - iprec ;
 111                 d_paths[idpath*(r-1)+p-1].y = 1 ;
 112           }
 113           iprec = ic ;
 114         }
 115         idpath++ ;
 116   }
 117   // Q1 sup
 118   for (int a=2 ; a>0 ; a--){         // les 2 angles 60 et 75
 119         for (int p=0 ; p< r ; p++){      // les r points
 120           jc = floor(tangente[a]*p + offset) ;
 121           if ( p > 0 ){
 122                 d_paths[idpath*(r-1)+p-1].x = -1 ;
 123                 d_paths[idpath*(r-1)+p-1].y = jc - jprec ;
 124           }
 125           jprec = jc ;
 126         }
 127         idpath++ ;
 128   }
 129
 130   // Q2
 131   basepath += 6 ;
 132   for (int a=0 ; a< 6 ; a++){         // les 6 angles 90,105,120,135,150,165
 133         for (int p=0 ; p<r-1 ; p++){      // les r points
 134           d_paths[idpath*(r-1)+p].x = -d_paths[(idpath - basepath)*(r-1)+p].y ;
 135           d_paths[idpath*(r-1)+p].y =  d_paths[(idpath - basepath)*(r-1)+p].x ;
 136           }
 137         idpath++ ;
 138   }
 139
 140   // Q3
 141   basepath += 6 ;
 142   for (int a=0 ; a< 6 ; a++){         // les 6 angles 180,195,210,225,240,255
 143         for (int p=0 ; p<r-1 ; p++){      // les r points
 144           d_paths[idpath*(r-1)+p].x = -d_paths[(idpath - basepath)*(r-1)+p].x ;
 145           d_paths[idpath*(r-1)+p].y = -d_paths[(idpath - basepath)*(r-1)+p].y ;
 146           }
 147         idpath++ ;
 148   }
 149
 150   // Q4
 151   basepath += 6 ;
 152   for (int a=0 ; a< 6 ; a++){         // les 6 angles 270,285,300,315,330,345
 153         for (int p=0 ; p<r-1 ; p++){      // les r points
 154           d_paths[idpath*(r-1)+p].x =  d_paths[(idpath - basepath)*(r-1)+p].y ;
 155           d_paths[idpath*(r-1)+p].y = -d_paths[(idpath - basepath)*(r-1)+p].x ;
 156           }
 157         idpath++ ;
 158   }
 159 }
 160
 161 /**
 162  *
 163  * \brief calcule l'estimation initiale
 164  * \author zulu - AND
 165  *
 166  * \param[in] L         Largeur de l'image
 167  * \param[in] H         Hauteur de l'image
 168  * \param[in] r         coté de la fenêtre de moyenneage
 169  *
 170  * \param[out] d_estim  Image estimee 0
 171  *
 172  * Version texture : l'img originale est supposée en texture.
 173  * L'estimation réalisée correspond a un moyenneur de 'rayon' r
 174  * Execution sur des blocs de threads 2D et une grille 2D
 175  * selon les dimensions de l'image.
 176  *
 177  */
 178 __global__ void kernel_neutre_img2estim(unsigned int *d_estim, unsigned int L, unsigned int H){
 179
 180   // coordonnes du point dans l'image
 181   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 182   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 183   unsigned int pos = i*L +j ;
 184
 185   d_estim[ pos ] = tex2D(tex_img_in, j, i) ;
 186
 187 }
 188
 189
 190 /**
 191  *
 192  * \brief calcule l'estimation initiale
 193  * \author zulu - AND
 194  *
 195  * \param[in] L         Largeur de l'image
 196  * \param[in] H         Hauteur de l'image
 197  * \param[in] r         coté de la fenêtre de moyenneage
 198  *
 199  * \param[out] d_estim  Image estimee 0
 200  *
 201  * Version texture : l'img originale est supposée en texture.
 202  * L'estimation réalisée correspond a un moyenneur de 'rayon' r
 203  * Execution sur des blocs de threads 2D et une grille 2D
 204  * selon les dimensions de l'image.
 205  *
 206  */
 207 __global__ void kernel_init_estim_from_img_in(unsigned int * d_estim, unsigned int L, unsigned int H, unsigned int r){
 208   // coordonnes du point dans l'image
 209   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 210   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 211   unsigned int pos = i*L +j ;
 212   unsigned int ic , jc, ng ;
 213
 214   if( (i>r)&&(i<H-r)&&(j>r)&&(j<L-r) ){
 215         ng = 0 ;
 216         for (ic = i - r ; ic <= i + r ; ic++){
 217           for (jc = j - r ; jc <= j + r ; jc++){
 218                 ng += tex2D(tex_img_in, jc, ic ) ;
 219           }
 220         }
 221         d_estim[ pos ] = ng/((2*r+1)*(2*r+1)) ;
 222   }
 223   /*
 224         // pour les bords  : pas de traitement
 225   else
 226   d_estim[ pos ] = tex2D(tex_img_in, j, i) ;
 227   */
 228 }
 229
 230
 231
 232 /**
 233  *
 234  * \brief calcule l'estimation initiale
 235  * \author zulu - AND
 236  *
 237  * \param[in] d_data    image originale
 238  * \param[in] L         Largeur de l'image
 239  * \param[in] H         Hauteur de l'image
 240  * \param[in] r         coté de la fenêtre de moyenneage
 241  *
 242  * \param[out] d_estim  Image estimee 0
 243  *
 244  * Version global mem : l'img originale est en mémoire globale, passée en param.
 245  * L'estimation réalisée correspond a un moyenneur de 'rayon' r
 246  * Execution sur des blocs de threads 2D et une grille 2D
 247  * selon les dimensions de l'image.
 248  * Moins rapide que les 2 autres solutions.
 249  *
 250  */
 251 __global__ void kernel_init_estim_from_img_in_global_mem(unsigned int * d_data, unsigned int * d_estim,
 252                                                                                                                  unsigned int L, unsigned int H, unsigned int r){
 253   // coordonnes du point dans l'image
 254   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 255   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 256   unsigned int pos = i*L +j ;
 257   unsigned int ic , jc, ng ;
 258
 259   if( (i>r)&&(i<H-r)&&(j>r)&&(j<L-r) ){
 260         ng = 0 ;
 261         for (ic = i - r ; ic <= i + r ; ic++){
 262           for (jc = j - r ; jc <= j + r ; jc++){
 263                 ng += d_data[ ic*L + jc ] ;
 264           }
 265         }
 266         d_estim[ pos ] = ng/((2*r+1)*(2*r+1)) ;
 267   }
 268 }
 269
 270
 271
 272 /**
 273  *
 274  * \brief calcule les niveaux de gris de l'estimation n+1
 275  * \author zulu - AND
 276  *
 277  * \param[in] L         Largeur de l'image
 278  * \param[in] H         Hauteur de l'image
 279  * \param[in] p         poids du terme lniv
 280  *
 281  * \param[out] d_estim  Image estimee n+1
 282  *
 283  * Version mixte : l'img originale est supposee en texture,
 284  * l'img lniv en mémoire globale, passée en param.
 285  * Cela évite la copie en texture de l'img lniv à chaque itération.
 286  * Execution sur des blocs de threads 2D et une grille 2D
 287  * selon les dimensions de l'image.
 288  * Moins rapide que 'texture' mais plus rapide que 'globalmem'
 289  *
 290  */
 291 __global__ void kernel_estim_next_step_hybrid(unsigned int * d_estim, unsigned int * d_lniv, unsigned int L, unsigned int H, unsigned int p){
 292   // coordonnes du point dans l'image
 293   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 294   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 295   unsigned int pos = i*L +j ;
 296
 297   d_estim[ pos ] = ( tex2D(tex_img_in, j, i) + p*d_lniv[ pos ] )/(1+p) ;
 298
 299 }
 300
 301 /**
 302  *
 303  * \brief calcule les niveaux de gris de l'estimation n+1
 304  * \author zulu - AND
 305  *
 306  * \param[in] L         Largeur de l'image
 307  * \param[in] H         Hauteur de l'image
 308  * \param[in] p         poids du terme lniv
 309  *
 310  * \param[out] d_estim  Image estimee n+1
 311  *
 312  * Version texture : Les donnees (img originale, img lniv) sont supposees en textures.
 313  * Execution sur des blocs de threads 2D et une grille 2D
 314  * selon les dimensions de l'image.
 315  * Plus rapide que les 2 autres solutions
 316  *
 317  */
 318 __global__ void kernel_estim_next_step_texture(unsigned int * d_estim, unsigned int L, unsigned int H, unsigned int p){
 319   // coordonnes du point dans l'image
 320   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 321   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 322   unsigned int pos = i*L +j ;
 323
 324   d_estim[ pos ] = ( tex2D(tex_img_in, j, i) + p*tex2D(tex_img_lniv, j, i ))/(1+p) ;
 325
 326 }
 327
 328 /**
 329  *
 330  * \brief calcule les niveaux de gris de l'estimation n+1
 331  * \author zulu - AND
 332  *
 333  * \param[in] d_lniv    Image des lniv n
 334  * \param[in] d_data    Image originale
 335  * \param[in] L         Largeur de l'image
 336  * \param[in] H         Hauteur de l'image
 337  * \param[in] p         poids du terme lniv
 338  *
 339  * \param[out] d_estim  Image estimee n+1
 340  *
 341  * Version mémoire globale : les données sont passées en params.
 342  * Execution sur des blocs de threads 2D et une grille 2D
 343  * selon les dimensions de l'image.
 344  *
 345  */
 346 __global__ void kernel_estim_next_step_global_mem(unsigned int * d_estim, unsigned int * d_lniv, unsigned int * d_data,
 347                                                                                                   unsigned int L, unsigned int H, unsigned int p){
 348   // coordonnes du point dans l'image
 349   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 350   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 351   unsigned int pos = i*L +j ;
 352
 353   d_estim[ pos ] = ( d_data[ pos ] + p*d_lniv[ pos ])/(1+p) ;
 354
 355 }
 356
 357 /**
 358  *
 359  * \brief determine les lniv en chaque point de l'image
 360  * \author zulu - AND
 361  *
 362  * \param[in] L         Largeur de l'image
 363  * \param[in] H         Hauteur de l'image
 364  * \param[in] r         longueur des segments
 365  *
 366  * \param[out] img_out  image des lniv
 367  *
 368  * Execution sur des blocs de threads 2D et une grille 2D
 369  * selon les dimensions de l'image.
 370  * L'image d'entrée doit être au préalable en mémoire texture pointée par "tex_img_estim".
 371  * Les matrices des chemins sont, elles, pointées par "tex_paths"
 372  * Cette version ne fournit pas les indices des chemins pour les tracé éventuel des lniv.
 373  */
 374 __global__ void kernel_levelines_texture(unsigned int * img_out, unsigned int L, unsigned int H, unsigned int r)
 375 {
 376   // coordonnes du point dans l'image
 377   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 378   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 379
 380   // nb de points par chemin
 381   int lpath =  r ;
 382   unsigned int ic, jc, zc, z ;
 383   int idpath, idpix ;
 384   unsigned int mse_min, mse_cur, val ;
 385   uint2 mse ;
 386
 387
 388   if((i>=lpath)&&(i<=H-lpath)&&(j>=lpath)&&(j<=L-lpath)){
 389         z = tex2D(tex_img_estim, j, i) ;
 390         for( idpath=0; idpath < PSIZE_I ; idpath++) {
 391           ic = i ;
 392           jc = j ;
 393           mse.x = z ;
 394           mse.y = z*z ;
 395           for( idpix=0; idpix < lpath-1 ; idpix++ ) {
 396                 ic += tex2D(tex_paths, idpix, idpath).x ;
 397                 jc += tex2D(tex_paths, idpix, idpath).y ;
 398                 zc = tex2D(tex_img_estim, jc, ic) ;
 399                 mse.x += zc ;
 400                 mse.y += zc*zc ;
 401           }
 402           // critere de selection du chemin ( SUM_(X2) - SUM_(X)2 / lpath )
 403           // TODO cherchera  ameliorer pour vitesse
 404           mse_cur = ( mse.y - ( mse.x / lpath ) * mse.x ) ;
 405           if ( (idpath == 0) || (mse_cur < mse_min) ) {
 406                 mse_min = mse_cur ;
 407                 val = mse.x ;
 408           }
 409         }
 410         img_out[ i*L + j ] = val / lpath ;
 411   }
 412 }
 413
 414
 415 /**
 416  *
 417  * \brief determine les lniv en chaque point de l'image
 418  * \author zulu - AND
 419  *
 420  * \param[in] L         Largeur de l'image
 421  * \param[in] H         Hauteur de l'image
 422  * \param[in] r         longueur des segments
 423  *
 424  * \param[out] img_out  image des lniv
 425  *
 426  * Execution sur des blocs de threads 2D et une grille 2D
 427  * selon les dimensions de l'image.
 428  * L'image d'entrée doit être au préalable en mémoire texture pointée par "tex_img_estim".
 429  * Les matrices des chemins sont, elles, pointées par "tex_paths"
 430  * Cette version ne fournit pas les indices des chemins pour les tracé éventuel des lniv.
 431  * Cette version tente d'utiliser la shared memory pour compenser la baisse de perf due aux chemins
 432  * paramétrables non constants.
 433  */
 434
 435 __global__ void kernel_levelines_texture_smem(unsigned int * img_out, unsigned int L, unsigned int H, unsigned int r)
 436 {
 437   // coordonnées du point dans le bloc
 438   unsigned int iib = threadIdx.x ;
 439   unsigned int jib = threadIdx.y ;
 440   // coordonnes du point dans l'image
 441   unsigned int i = blockIdx.x*blockDim.x + iib ;
 442   unsigned int j = blockIdx.y*blockDim.y + jib ;
 443
 444   // nb de points par chemin
 445   int lpath =  r ;
 446   int ic, jc ;
 447   int idpath, idpix ;
 448   unsigned int val, mse_cur, mse_min, z, zc ;
 449   uint2 mse_data ;
 450
 451   //__shared__ unsigned int val_img[16*16] ;
 452
 453   //val_img[jib*16+iib] = tex2D(tex_img_estim, j, i) ;
 454
 455   if((i>=lpath)&&(i<=H-lpath)&&(j>=lpath)&&(j<=L-lpath)){
 456         z = tex2D(tex_img_estim, j, i) ;
 457         for( idpath=0; idpath < PSIZE_I ; idpath++) {
 458           ic = i ;
 459           jc = j ;
 460           mse_data.x = z ;
 461           mse_data.y = z*z ;
 462           mse_min = mse_data.y - mse_data.x/lpath*mse_data.y ;
 463           for( idpix=0; idpix < lpath-1 ; idpix++ ) {
 464                 ic += tex2D(tex_paths, idpix, idpath).x ;
 465                 jc += tex2D(tex_paths, idpix, idpath).y ;
 466                 zc = tex2D(tex_img_estim, jc, ic) ;
 467                 mse_data.x += zc ;
 468                 mse_data.y += zc*zc ;
 469           }
 470           // critere de selection du chemin ( SUM_(X2) - SUM_(X)2 / lpath )
 471           // TODO cherchera  ameliorer pour vitesse
 472           mse_cur = ( mse_data.y - ( mse_data.x / lpath ) * mse_data.x ) ;
 473           if ( mse_cur < mse_min ){
 474                 mse_min = mse_cur ;
 475                 val = mse_data.x ;
 476           }
 477         }
 478         img_out[ i*L + j ] = val / lpath ;
 479   }
 480 }
 481
 482 /**
 483  *
 484  * \brief trace les segments sur un maillage carré
 485  * \author zulu - AND
 486  *
 487  * \param[in] img_in    image d'entree
 488  * \param[in] dir       tableaux des directions
 489  * \param[in] L         Largeur de l'image
 490  * \param[in] H         Hauteur de l'image
 491  * \param[in] pas       coté du maillage
 492  * \param[in] ng        niveau de gris des segments
 493  * \param[in] r         longueur des segments
 494  *
 495  * \param[out] img_out  image + les segments superposés
 496  *
 497  * Kernel trivial. Ne trace rien sur les bords.
 498  * execution sur des blocs de threads 2D et une grille 2D
 499  * selon les dimensions de l'image
 500  */
 501 __global__ void kernel_trace_levelines(unsigned int * img_in, unsigned int * dir, unsigned int * img_out,
 502                                                                            unsigned int L, unsigned int H, unsigned int pas, unsigned int ng,
 503                                                                            unsigned int r ){
 504   // coordonnes du point dans l'image
 505   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 506   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;;
 507
 508   // nb de points par chemin
 509   int lpath =  r ;
 510   unsigned int ic, jc, idpix ;
 511   unsigned int idpath  ;
 512
 513   img_out[ i*L + j ] = img_in[ i*L + j ] ;
 514
 515   if ( !(i%pas+j%pas)&&(i>lpath)&&(i<H-lpath)&&(j>lpath)&&(j<L-lpath) ){
 516         ic = i ;
 517         jc = j ;
 518         idpath = dir[ic*L+jc] ;
 519         img_out[ ic*L+jc ] = ng ;
 520         for ( idpix=0 ; idpix < lpath-1 ; idpix++ ){
 521           ic += tex2D(tex_paths, idpix, idpath).x ; // pathDi[idpath][idpix] ;
 522           jc += tex2D(tex_paths, idpix, idpath).y ; // pathDj[idpath][idpix] ;
 523           img_out[ ic*L + jc ] = ng ;
 524           }
 525   }
 526
 527 }