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<ushort, 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( ushort * 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   char MSQ, LSQ ;
 105
 106   // Q1 inf
 107   for (int a=0 ; a< 4 ; a++){        // les 4 angles 0,15,30 et 45
 108         for (int p=0 ; p< r ; p++){      // les r points
 109           ic = r-1 - floor(tangente[a]*p + offset) ;
 110           if ( p > 0 ){
 111                 MSQ = ic - iprec ;
 112                 LSQ = 1 ;
 113                 //d_paths[idpath*(r-1)+p-1].x = ic - iprec ;
 114                 //d_paths[idpath*(r-1)+p-1].y = 1 ;
 115                 d_paths[idpath*(r-1)+p-1] = ((short)MSQ << 8) | LSQ ;
 116           }
 117           iprec = ic ;
 118         }
 119         idpath++ ;
 120   }
 121   // Q1 sup
 122   for (int a=2 ; a>0 ; a--){         // les 2 angles 60 et 75
 123         for (int p=0 ; p< r ; p++){      // les r points
 124           jc = floor(tangente[a]*p + offset) ;
 125           if ( p > 0 ){
 126                 MSQ = -1 ;
 127                 LSQ = jc - jprec ;
 128                 d_paths[idpath*(r-1)+p-1] = ((short)MSQ << 8) | LSQ ;
 129                 //d_paths[idpath*(r-1)+p-1].x = -1 ;
 130                 //d_paths[idpath*(r-1)+p-1].y = jc - jprec ;
 131           }
 132           jprec = jc ;
 133         }
 134         idpath++ ;
 135   }
 136
 137   // Q2
 138   basepath += 6 ;
 139   for (int a=0 ; a< 6 ; a++){         // les 6 angles 90,105,120,135,150,165
 140         for (int p=0 ; p<r-1 ; p++){      // les r points
 141           MSQ = - ( d_paths[(idpath - basepath)*(r-1)+p] & 0x00FF ) ;
 142           LSQ = ( d_paths[(idpath - basepath)*(r-1)+p] >> 8 ) ;
 143           d_paths[idpath*(r-1)+p-1] = ((short)MSQ << 8) | LSQ ;
 144           //d_paths[idpath*(r-1)+p].x = -d_paths[(idpath - basepath)*(r-1)+p].y ;
 145           //d_paths[idpath*(r-1)+p].y =  d_paths[(idpath - basepath)*(r-1)+p].x ;
 146           }
 147         idpath++ ;
 148   }
 149
 150   // Q3
 151   basepath += 6 ;
 152   for (int a=0 ; a< 6 ; a++){         // les 6 angles 180,195,210,225,240,255
 153         for (int p=0 ; p<r-1 ; p++){      // les r points
 154           MSQ = - ( d_paths[(idpath - basepath)*(r-1)+p] >> 8 ) ;
 155           LSQ = - ( d_paths[(idpath - basepath)*(r-1)+p] & 0x00FF ) ;
 156           d_paths[idpath*(r-1)+p-1] = ((short)MSQ << 8) | LSQ ;
 157           //d_paths[idpath*(r-1)+p].x = -d_paths[(idpath - basepath)*(r-1)+p].x ;
 158           //d_paths[idpath*(r-1)+p].y = -d_paths[(idpath - basepath)*(r-1)+p].y ;
 159           }
 160         idpath++ ;
 161   }
 162
 163   // Q4
 164   basepath += 6 ;
 165   for (int a=0 ; a< 6 ; a++){         // les 6 angles 270,285,300,315,330,345
 166         for (int p=0 ; p<r-1 ; p++){      // les r points
 167           MSQ = d_paths[(idpath - basepath)*(r-1)+p] & 0x00FF ;
 168           LSQ = - ( d_paths[(idpath - basepath)*(r-1)+p] >> 8 ) ;
 169           d_paths[idpath*(r-1)+p-1] = ((short)MSQ << 8) | LSQ ;
 170           //d_paths[idpath*(r-1)+p].x =  d_paths[(idpath - basepath)*(r-1)+p].y ;
 171           //d_paths[idpath*(r-1)+p].y = -d_paths[(idpath - basepath)*(r-1)+p].x ;
 172           }
 173         idpath++ ;
 174   }
 175
 176 }
 177
 178 /**
 179  *
 180  * \brief calcule l'estimation initiale
 181  * \author zulu - AND
 182  *
 183  * \param[in] L         Largeur de l'image
 184  * \param[in] H         Hauteur de l'image
 185  * \param[in] r         coté de la fenêtre de moyenneage
 186  *
 187  * \param[out] d_estim  Image estimee 0
 188  *
 189  * Version texture : l'img originale est supposée en texture.
 190  * L'estimation réalisée correspond a un moyenneur de 'rayon' r
 191  * Execution sur des blocs de threads 2D et une grille 2D
 192  * selon les dimensions de l'image.
 193  *
 194  */
 195 __global__ void kernel_neutre_img2estim(unsigned int *d_estim, unsigned int L, unsigned int H){
 196
 197   // coordonnes du point dans l'image
 198   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 199   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 200   unsigned int pos = i*L +j ;
 201
 202   d_estim[ pos ] = tex2D(tex_img_in, j, i) ;
 203
 204 }
 205
 206
 207 /**
 208  *
 209  * \brief calcule l'estimation initiale
 210  * \author zulu - AND
 211  *
 212  * \param[in] L         Largeur de l'image
 213  * \param[in] H         Hauteur de l'image
 214  * \param[in] r         coté de la fenêtre de moyenneage
 215  *
 216  * \param[out] d_estim  Image estimee 0
 217  *
 218  * Version texture : l'img originale est supposée en texture.
 219  * L'estimation réalisée correspond a un moyenneur de 'rayon' r
 220  * Execution sur des blocs de threads 2D et une grille 2D
 221  * selon les dimensions de l'image.
 222  *
 223  */
 224 __global__ void kernel_init_estim_from_img_in(unsigned int * d_estim, unsigned int L, unsigned int H, unsigned int r){
 225   // coordonnes du point dans l'image
 226   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 227   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 228   unsigned int pos = i*L +j ;
 229   unsigned int ic , jc, ng ;
 230
 231   if( (i>r)&&(i<H-r)&&(j>r)&&(j<L-r) ){
 232         ng = 0 ;
 233         for (ic = i - r ; ic <= i + r ; ic++){
 234           for (jc = j - r ; jc <= j + r ; jc++){
 235                 ng += tex2D(tex_img_in, jc, ic ) ;
 236           }
 237         }
 238         d_estim[ pos ] = ng/((2*r+1)*(2*r+1)) ;
 239   }
 240   /*
 241         // pour les bords  : pas de traitement
 242   else
 243   d_estim[ pos ] = tex2D(tex_img_in, j, i) ;
 244   */
 245 }
 246
 247
 248
 249 /**
 250  *
 251  * \brief calcule l'estimation initiale
 252  * \author zulu - AND
 253  *
 254  * \param[in] d_data    image originale
 255  * \param[in] L         Largeur de l'image
 256  * \param[in] H         Hauteur de l'image
 257  * \param[in] r         coté de la fenêtre de moyenneage
 258  *
 259  * \param[out] d_estim  Image estimee 0
 260  *
 261  * Version global mem : l'img originale est en mémoire globale, passée en param.
 262  * L'estimation réalisée correspond a un moyenneur de 'rayon' r
 263  * Execution sur des blocs de threads 2D et une grille 2D
 264  * selon les dimensions de l'image.
 265  * Moins rapide que les 2 autres solutions.
 266  *
 267  */
 268 __global__ void kernel_init_estim_from_img_in_global_mem(unsigned int * d_data, unsigned int * d_estim,
 269                                                                                                                  unsigned int L, unsigned int H, unsigned int r){
 270   // coordonnes du point dans l'image
 271   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 272   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 273   unsigned int pos = i*L +j ;
 274   unsigned int ic , jc, ng ;
 275
 276   if( (i>r)&&(i<H-r)&&(j>r)&&(j<L-r) ){
 277         ng = 0 ;
 278         for (ic = i - r ; ic <= i + r ; ic++){
 279           for (jc = j - r ; jc <= j + r ; jc++){
 280                 ng += d_data[ ic*L + jc ] ;
 281           }
 282         }
 283         d_estim[ pos ] = ng/((2*r+1)*(2*r+1)) ;
 284   }
 285 }
 286
 287
 288
 289 /**
 290  *
 291  * \brief calcule les niveaux de gris de l'estimation n+1
 292  * \author zulu - AND
 293  *
 294  * \param[in] L         Largeur de l'image
 295  * \param[in] H         Hauteur de l'image
 296  * \param[in] p         poids du terme lniv
 297  *
 298  * \param[out] d_estim  Image estimee n+1
 299  *
 300  * Version mixte : l'img originale est supposee en texture,
 301  * l'img lniv en mémoire globale, passée en param.
 302  * Cela évite la copie en texture de l'img lniv à chaque itération.
 303  * Execution sur des blocs de threads 2D et une grille 2D
 304  * selon les dimensions de l'image.
 305  * Moins rapide que 'texture' mais plus rapide que 'globalmem'
 306  *
 307  */
 308 __global__ void kernel_estim_next_step_hybrid(unsigned int * d_estim, unsigned int * d_lniv, unsigned int L, unsigned int H, unsigned int p){
 309   // coordonnes du point dans l'image
 310   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 311   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 312   unsigned int pos = i*L +j ;
 313
 314   d_estim[ pos ] = ( tex2D(tex_img_in, j, i) + p*d_lniv[ pos ] )/(1+p) ;
 315
 316 }
 317
 318 /**
 319  *
 320  * \brief calcule les niveaux de gris de l'estimation n+1
 321  * \author zulu - AND
 322  *
 323  * \param[in] L         Largeur de l'image
 324  * \param[in] H         Hauteur de l'image
 325  * \param[in] p         poids du terme lniv
 326  *
 327  * \param[out] d_estim  Image estimee n+1
 328  *
 329  * Version texture : Les donnees (img originale, img lniv) sont supposees en textures.
 330  * Execution sur des blocs de threads 2D et une grille 2D
 331  * selon les dimensions de l'image.
 332  * Plus rapide que les 2 autres solutions
 333  *
 334  */
 335 __global__ void kernel_estim_next_step_texture(unsigned int * d_estim, unsigned int L, unsigned int H, unsigned int p){
 336   // coordonnes du point dans l'image
 337   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 338   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 339   unsigned int pos = i*L +j ;
 340
 341   d_estim[ pos ] = ( tex2D(tex_img_in, j, i) + p*tex2D(tex_img_lniv, j, i ))/(1+p) ;
 342
 343 }
 344
 345 /**
 346  *
 347  * \brief calcule les niveaux de gris de l'estimation n+1
 348  * \author zulu - AND
 349  *
 350  * \param[in] d_lniv    Image des lniv n
 351  * \param[in] d_data    Image originale
 352  * \param[in] L         Largeur de l'image
 353  * \param[in] H         Hauteur de l'image
 354  * \param[in] p         poids du terme lniv
 355  *
 356  * \param[out] d_estim  Image estimee n+1
 357  *
 358  * Version mémoire globale : les données sont passées en params.
 359  * Execution sur des blocs de threads 2D et une grille 2D
 360  * selon les dimensions de l'image.
 361  *
 362  */
 363 __global__ void kernel_estim_next_step_global_mem(unsigned int * d_estim, unsigned int * d_lniv, unsigned int * d_data,
 364                                                                                                   unsigned int L, unsigned int H, unsigned int p){
 365   // coordonnes du point dans l'image
 366   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 367   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 368   unsigned int pos = i*L +j ;
 369
 370   d_estim[ pos ] = ( d_data[ pos ] + p*d_lniv[ pos ])/(1+p) ;
 371
 372 }
 373
 374 /**
 375  *
 376  * \brief determine les lniv en chaque point de l'image
 377  * \author zulu - AND
 378  *
 379  * \param[in] L         Largeur de l'image
 380  * \param[in] H         Hauteur de l'image
 381  * \param[in] r         longueur des segments
 382  *
 383  * \param[out] img_out  image des lniv
 384  *
 385  * Execution sur des blocs de threads 2D et une grille 2D
 386  * selon les dimensions de l'image.
 387  * L'image d'entrée doit être au préalable en mémoire texture pointée par "tex_img_estim".
 388  * Les matrices des chemins sont, elles, préalablement chargées en SHMEM depuis la texture"
 389  * Cette version ne fournit pas les indices des chemins pour les tracé éventuel des lniv.
 390  */
 391 __global__ void kernel_levelines_texture(unsigned int * img_out, unsigned int L, unsigned int H, unsigned int r)
 392 {
 393   // coordonnees du point dans le bloc
 394   unsigned int iib = threadIdx.x ;
 395   unsigned int jib = threadIdx.y ;
 396   // coordonnees du point dans l'image
 397   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 398   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 399
 400   // nb de points par chemin
 401   int lpath =  r ;
 402   unsigned int ic, jc, zc, z ;
 403   int idpath, idpix ;
 404   unsigned int mse_min, mse_cur, val ;
 405   uint2 mse ;
 406   short texVal ;
 407
 408
 409   extern __shared__ short shPath[] ;
 410
 411   unsigned int absPos = jib*8 + iib ;
 412   if ( absPos < PSIZE_I ){
 413         for ( idpix = 0; idpix < lpath-1; idpix++){
 414           shPath[ idpix*24 + absPos ] = tex2D(tex_paths, idpix, absPos) ;
 415         }
 416         syncthreads() ;
 417   }
 418
 419   if((i>=lpath)&&(i<=H-lpath)&&(j>=lpath)&&(j<=L-lpath)){
 420         z = tex2D(tex_img_estim, j, i) ;
 421         for( idpath=0; idpath < PSIZE_I ; idpath++) {
 422           ic = i ;
 423           jc = j ;
 424           mse.x = z ;
 425           mse.y = z*z ;
 426           for( idpix=0; idpix < lpath-1 ; idpix++ ) {
 427                 texVal = shPath[ idpix*24 + idpath ] ;
 428                 ic += (char)(texVal>>8) ;
 429                 jc += (char)(texVal) ;
 430                 zc = tex2D(tex_img_estim, jc, ic) ;
 431                 mse.x += zc ;
 432                 mse.y += zc*zc ;
 433           }
 434           // critere de selection du chemin ( SUM_(X2) - SUM_(X)2 / lpath )
 435           // TODO cherchera  ameliorer pour vitesse
 436           mse_cur = ( mse.y - ( mse.x / lpath ) * mse.x ) ;
 437           if ( (idpath == 0) || (mse_cur < mse_min) ) {
 438                 mse_min = mse_cur ;
 439                 val = mse.x ;
 440           }
 441         }
 442         img_out[ i*L + j ] = val / lpath ;
 443   }
 444 }
 445
 446
 447 /**
 448  *
 449  * \brief trace les segments sur un maillage carré
 450  * \author zulu - AND
 451  *
 452  * \param[in] img_in    image d'entree
 453  * \param[in] dir       tableaux des directions
 454  * \param[in] L         Largeur de l'image
 455  * \param[in] H         Hauteur de l'image
 456  * \param[in] pas       coté du maillage
 457  * \param[in] ng        niveau de gris des segments
 458  * \param[in] r         longueur des segments
 459  *
 460  * \param[out] img_out  image + les segments superposés
 461  *
 462  * Kernel trivial. Ne trace rien sur les bords.
 463  * execution sur des blocs de threads 2D et une grille 2D
 464  * selon les dimensions de l'image
 465  */
 466 /*
 467 __global__ void kernel_trace_levelines(unsigned int * img_in, unsigned int * dir, unsigned int * img_out,
 468                                                                            unsigned int L, unsigned int H, unsigned int pas, unsigned int ng,
 469                                                                            unsigned int r ){
 470   // coordonnes du point dans l'image
 471   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 472   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;;
 473
 474   // nb de points par chemin
 475   int lpath =  r ;
 476   unsigned int ic, jc, idpix ;
 477   unsigned int idpath  ;
 478
 479   img_out[ i*L + j ] = img_in[ i*L + j ] ;
 480
 481   if ( !(i%pas+j%pas)&&(i>lpath)&&(i<H-lpath)&&(j>lpath)&&(j<L-lpath) ){
 482         ic = i ;
 483         jc = j ;
 484         idpath = dir[ic*L+jc] ;
 485         img_out[ ic*L+jc ] = ng ;
 486         for ( idpix=0 ; idpix < lpath-1 ; idpix++ ){
 487           ic += tex2D(tex_paths, idpix, idpath).x ; // pathDi[idpath][idpix] ;
 488           jc += tex2D(tex_paths, idpix, idpath).y ; // pathDj[idpath][idpix] ;
 489           img_out[ ic*L + jc ] = ng ;
 490           }
 491   }
 492
 493 }
 494 */