src/lib_kernel_snake_2_gpu.cu

   1
   2
   3 __global__ void genere_snake_rectangle_4nodes_gpu(snake_node_gpu * d_snake, int dist_bords, int i_dim, int j_dim){
   4   if (threadIdx.x == 0){
   5         int n = 0;
   6         /* n0 */
   7         d_snake[n].posi = dist_bords ;
   8         d_snake[n].posj = dist_bords ;
   9         n++ ;
  10         /* n1 */
  11         d_snake[n].posi = i_dim - dist_bords ;
  12         d_snake[n].posj = dist_bords ;
  13         n++ ;
  14         /* n2 */
  15         d_snake[n].posi = i_dim - dist_bords ;
  16         d_snake[n].posj = j_dim - dist_bords ;
  17         n++ ;
  18         /* n3 */
  19         d_snake[n].posi = dist_bords ;
  20         d_snake[n].posj = j_dim - dist_bords ;
  21
  22         for (int i=0; i<4; i++)
  23           {
  24                 d_snake[i].freeman_in = 0;
  25                 d_snake[i].freeman_out = 0;
  26                 d_snake[i].centre_i = 0;
  27                 d_snake[i].centre_j = 0;
  28                 d_snake[i].last_move = 0;
  29                 d_snake[i].nb_pixels = 123;
  30                 d_snake[i].code_segment = 0;
  31
  32           }
  33   }
  34 }
  35
  36 __global__ void genere_diagos_rectangle(uint4 * d_diagos, int h, int l, int q){
  37   int inci = h/q;
  38   int incj = l/q;
  39   int iM,jM, iN, jN ;
  40   int idxDiago = 0 ;
  41         // boucle double pour les positions du point NO de la diagonale
  42         for ( iM = 0; iM < q-1; iM++){
  43           for (jM = 0 ; jM < q-1 ; jM++){
  44                 //boucle double pour les positions du point SE de la diagonale
  45                 for (iN = iM+1; iN < q; iN++){
  46                   for (jN = jM+1; jN < q; jN++){
  47                         d_diagos[idxDiago].x = iM*inci;
  48                         d_diagos[idxDiago].y = jM*incj;
  49                         d_diagos[idxDiago].z = iN*inci;
  50                         d_diagos[idxDiago].w = jN*incj;
  51                         idxDiago++;
  52                   }
  53                 }
  54           }
  55         }
  56 }
  57
  58 __global__ void genere_snake_rectangle_Nnodes_gpu(snake_node_gpu * d_snake, int dist_bords, int i_dim, int j_dim){
  59   int nb_node_seg = 9 ;
  60   int limite = 64 ;
  61
  62   int i , h= i_dim-2*dist_bords, l= j_dim-2*dist_bords ;
  63   int inch = h/(nb_node_seg+1), incl= l/(nb_node_seg+1) ;
  64   if (threadIdx.x == 0){
  65         int n = 0;
  66         /* n0 */
  67         d_snake[n].posi = dist_bords ;
  68         d_snake[n].posj = dist_bords ;
  69         n++ ;
  70         /*entre sommet 0 et 1*/
  71         i = 0 ;
  72         while (i < nb_node_seg)
  73           {
  74                 if ( (d_snake[n-1].posi + inch)-(i_dim - dist_bords) > limite )
  75                   d_snake[n].posi = d_snake[n-1].posi + inch ;
  76                 else
  77                   d_snake[n].posi = d_snake[n-1].posi + inch/2 ;
  78                 d_snake[n].posj = dist_bords ;
  79                 d_snake[n-1].nb_pixels =  d_snake[n].posi -  d_snake[n-1].posi ;
  80                 n++ ; i++ ;
  81           }
  82         /* n1 */
  83         d_snake[n].posi = i_dim - dist_bords ;
  84         d_snake[n].posj = dist_bords ;
  85         d_snake[n-1].nb_pixels =  d_snake[n].posi -  d_snake[n-1].posi ;
  86         n++ ;
  87         /*entre S1 et S2*/
  88         i = 0 ;
  89         while (i< nb_node_seg)
  90           {
  91                 if ( (j_dim - dist_bords) - (d_snake[n-1].posj + incl) > limite )
  92                   d_snake[n].posj = d_snake[n-1].posj + incl ;
  93                 else
  94                   d_snake[n].posj = d_snake[n-1].posj + incl/2 ;
  95                 d_snake[n].posi = i_dim - dist_bords ;
  96                 d_snake[n-1].nb_pixels =  d_snake[n].posj -  d_snake[n-1].posj ;
  97                 n++ ; i++ ;
  98           }
  99         /* n2 */
 100         d_snake[n].posi = i_dim - dist_bords ;
 101         d_snake[n].posj = j_dim - dist_bords ;
 102         d_snake[n-1].nb_pixels =  d_snake[n].posj -  d_snake[n-1].posj ;
 103         n++ ;
 104         /*entre S2 et S3*/
 105         i = 0 ;
 106         while (i< nb_node_seg)
 107           {
 108                 if ( (d_snake[n-1].posi - inch) - dist_bords > limite )
 109                   d_snake[n].posi = d_snake[n-1].posi - inch ;
 110                 else
 111                   d_snake[n].posi = d_snake[n-1].posi - inch/2 ;
 112                 d_snake[n].posj = j_dim - dist_bords ;
 113                 d_snake[n-1].nb_pixels =  d_snake[n-1].posi -  d_snake[n].posi ;
 114                 n++ ; i++ ;
 115           }
 116         /* n3 */
 117         d_snake[n].posi = dist_bords ;
 118         d_snake[n].posj = j_dim - dist_bords ;
 119         d_snake[n-1].nb_pixels =  d_snake[n-1].posi -  d_snake[n].posi ;
 120         n++ ;
 121         /*entre S3 et S0*/
 122         i = 0 ;
 123         while (i< nb_node_seg)
 124           {
 125                 if ( (d_snake[n-1].posj - incl) - dist_bords > limite)
 126                   d_snake[n].posj = d_snake[n-1].posj - incl ;
 127                 else
 128                   d_snake[n].posj = d_snake[n-1].posj - incl/2 ;
 129                 d_snake[n].posi = dist_bords ;
 130                 d_snake[n-1].nb_pixels =  d_snake[n-1].posj -  d_snake[n].posj ;
 131                 n++ ; i++ ;
 132           }
 133         d_snake[n-1].nb_pixels =  d_snake[n-1].posj -  d_snake[0].posj ;
 134         for (i=0; i<n; i++)
 135           {
 136                 d_snake[i].freeman_in = 0;
 137                 d_snake[i].freeman_out = 0;
 138                 d_snake[i].centre_i = 0;
 139                 d_snake[i].centre_j = 0;
 140                 d_snake[i].last_move = 1;
 141                 d_snake[i].code_segment = 0;
 142
 143           }
 144   }
 145 }
 146
 147 __global__ void calcul_contribs_segments_snake(snake_node_gpu * d_snake, int nb_nodes,
 148                                                                                            t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2,
 149                                                                                            int l, uint2 * liste_pix, t_sum_x2 * gsombloc, int * d_table_freeman)
 150 {
 151   // indices des elements
 152   int blockSize = blockDim.x ;
 153   int tib = threadIdx.x ;
 154   int nblocs_seg =  gridDim.x / nb_nodes ;
 155   int idx = blockDim.x*blockIdx.x + threadIdx.x ;
 156   int segment = blockIdx.x / nblocs_seg ;
 157   int tis = idx - segment*nblocs_seg*blockDim.x ;
 158
 159   //tab pour coordonnées pixels & contribs pixels de taille = (blockDim.x+offset(dec,dec2) )*(sizeof(t_sum_1+t_sum_x+t_sum_x2))
 160   extern __shared__ t_sum_1 scumuls_1[] ; // blockDim varie selon la longueur des segments => taille smem dynamique
 161   t_sum_x* scumuls_x = (t_sum_x*) &scumuls_1[CFI(blockDim.x)] ;
 162   t_sum_x2* scumuls_x2 = (t_sum_x2*) &scumuls_x[CFI(blockDim.x)] ;
 163
 164   //indices des noeuds
 165   uint x1, y1, x2, y2 ;
 166   int n1, n2 ;
 167
 168   n1 = segment ;
 169   n2 = segment +1 ;
 170   //gestion du bouclage du snake
 171   if (n2 >= nb_nodes) n2 = 0 ;
 172
 173   //affectation des differentes positions aux différents segments 'blocs de threads'
 174         x1 = d_snake[n1].posj ;
 175         y1 = d_snake[n1].posi ;
 176         x2 = d_snake[n2].posj ;
 177         y2 = d_snake[n2].posi ;
 178
 179   //params des deplacements
 180   int dx=x2-x1;
 181   int dy=y2-y1;
 182   uint abs_dx = ABS(dx);
 183   uint abs_dy = ABS(dy);
 184   uint nb_pix = abs_dy>abs_dx?(abs_dy+1):(abs_dx+1); // alternative -> lecture ds liste_points[]
 185   int incx=0, incy=0;
 186   uint2 p ;
 187   int xprec, xsuiv ;
 188
 189   //calcul liste des pixels du segment (x1,y1)-(x2,y2)
 190   if (dy > 0) incy=1; else incy=-1 ;
 191   if (dx > 0) incx=1; else incx=-1 ;
 192
 193   if (tis < nb_pix){
 194         if (abs_dy > abs_dx){
 195           //1 thread par ligne
 196           double k = (double)dx/dy ;
 197           p.x = y1 + incy*tis ;
 198           p.y =  x1 + floor((double)incy*k*tis+0.5) ;
 199           //enreg. coords. pixels en global mem pour freemans
 200
 201           if ((tis < 2)||(tis > nb_pix - 3)||(tis == nb_pix/2))
 202                 {
 203                   liste_pix[idx].x = p.x ;
 204                   liste_pix[idx].y = p.y ;
 205                   }
 206
 207         } else {
 208           //1 thread par colonne
 209           double k=(double)dy/dx ;
 210           p.x = y1 + floor((double)(incx*k*tis)+0.5) ;
 211           p.y = x1 + incx*tis ;
 212           if ( tis > 0 ){
 213                 xsuiv = y1 + floor((double)(incx*k*(tis+1))+0.5) ;
 214                 xprec = y1 + floor((double)(incx*k*(tis-1))+0.5) ;
 215           }
 216           //enreg. coords. pixels en global mem pour freeman
 217           //TODO
 218           //on peut calculer les freemans des segments
 219           //sans stocker l'ensemble des valeurs des pixels
 220           //juste avec les derivees aux extremites a calculer ici
 221
 222           if ((tis < 2)||(tis > nb_pix - 3)||(tis == nb_pix/2))
 223                 {
 224                   liste_pix[idx].x = p.x ;
 225                   liste_pix[idx].y = p.y ;
 226                   }
 227
 228         }
 229   }
 230   __syncthreads();
 231
 232   //calcul contribs individuelles des pixels
 233
 234   if ( (tis >0) && (tis < nb_pix-1)
 235            && ( (abs_dy <= abs_dx)
 236                         && ( (xprec > p.x) || (xsuiv > p.x))
 237                         || (abs_dy > abs_dx) ) )
 238         {
 239         int pos = p.x * l + p.y ;
 240         scumuls_1[ CFI(tib)] = 1+p.y;
 241         scumuls_x[ CFI(tib)] = cumul_x[ pos ] ;
 242         scumuls_x2[CFI(tib)] = cumul_x2[ pos ];
 243   } else {
 244         scumuls_1[ CFI(tib)] = 0;
 245         scumuls_x[ CFI(tib)] = 0;
 246         scumuls_x2[CFI(tib)] = 0;
 247   }
 248
 249    __syncthreads();
 250   //somme des contribs individuelles
 251   // unroll des  sommes partielles en smem
 252
 253   if (blockSize >= 512) {
 254         if (tib < 256) {
 255           scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib + 256) ];
 256           scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib + 256) ];
 257           scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib + 256) ];
 258         }
 259         __syncthreads();
 260   }
 261
 262   if (blockSize >= 256) {
 263         if (tib < 128) {
 264           scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib + 128) ];
 265           scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib + 128) ];
 266           scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib + 128) ];
 267         }
 268         __syncthreads();
 269   }
 270   if (blockSize >= 128) {
 271         if (tib <  64) {
 272           scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib +  64) ];
 273           scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib +  64) ];
 274           scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib +  64) ];
 275         }
 276         __syncthreads();
 277   }
 278
 279   //32 threads <==> 1 warp
 280   if (tib < 32)
 281         {
 282           {
 283                 scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib + 32) ];
 284                 scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib + 32) ];
 285                 scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib + 32) ];
 286           }
 287           {
 288                 scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib + 16) ];
 289                 scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib + 16) ];
 290                 scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib + 16) ];
 291           }
 292           {
 293                 scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib +  8) ];
 294                 scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib +  8) ];
 295                 scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib +  8) ];
 296           }
 297           scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib +  4) ];
 298           scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib +  4) ];
 299           scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib +  4) ];
 300
 301           scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib +  2) ];
 302           scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib +  2) ];
 303           scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib +  2) ];
 304
 305           scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib +  1) ];
 306           scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib +  1) ];
 307           scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib +  1) ];
 308         }
 309
 310   // resultat sommes partielles en gmem
 311   if (tib == 0) {
 312         gsombloc[ blockIdx.x ] = (t_sum_x2) scumuls_1[0];
 313         gsombloc[ blockIdx.x + gridDim.x ] = (t_sum_x2) scumuls_x[0];
 314         gsombloc[ blockIdx.x + 2*gridDim.x ] = (t_sum_x2) scumuls_x2[0];
 315   }
 316
 317   //calculs freemans, centre et code segment
 318   //1 uint4 par segment
 319
 320   int Di, Dj;
 321   if (tis == 0){
 322         Di = 1 + liste_pix[idx+1].x - liste_pix[idx].x ;
 323         Dj = 1 + liste_pix[idx+1].y - liste_pix[idx].y ;
 324         d_snake[segment].freeman_out = d_table_freeman[3*Di + Dj] ;
 325         //code seg
 326         if (dy > 0 ) d_snake[segment].code_segment = -1 ;
 327         if (dy < 0 ) d_snake[segment].code_segment = 1 ;
 328         if (dy == 0) d_snake[segment].code_segment = 0 ;
 329   }
 330
 331   if (tis == nb_pix-1){
 332         Di = 1 + liste_pix[idx].x - liste_pix[idx-1].x  ;
 333         Dj = 1 + liste_pix[idx].y - liste_pix[idx-1].y;
 334         d_snake[segment].freeman_in = d_table_freeman[3*Di + Dj] ;
 335   }
 336
 337   if (tis == (nb_pix/2)){
 338         d_snake[segment].centre_i = liste_pix[idx].x ;
 339         d_snake[segment].centre_j = liste_pix[idx].y ;
 340         }
 341 }
 342
 343 /*
 344   sommme des contribs par bloc -> contribs segment, pour le snake
 345
 346   execution sur : 1bloc / 1 thread par segment
 347  */
 348
 349 __global__ void somsom_snake(t_sum_x2 * somblocs, int nb_nodes, unsigned int nb_bl_seg, snake_node_gpu * d_snake){
 350
 351   t_sum_x2 sdata[3];
 352   unsigned int seg = blockIdx.x ;
 353
 354   //un thread par segment
 355   {
 356         sdata[0] = 0;
 357         sdata[1] = 0;
 358         sdata[2] = 0;
 359   }
 360
 361   for (int b=0; b < nb_bl_seg ; b++){
 362         sdata[0] += somblocs[seg*nb_bl_seg + b];
 363         sdata[1] += somblocs[(seg + nb_nodes)*nb_bl_seg + b];
 364         sdata[2] += somblocs[(seg + 2*nb_nodes)*nb_bl_seg + b];
 365   }
 366
 367   //totaux en gmem
 368   {
 369         d_snake[seg].sum_1 = sdata[0];
 370         d_snake[seg].sum_x = sdata[1];
 371         d_snake[seg].sum_x2 = sdata[2];
 372   }
 373 }
 374
 375 __device__ double codage_gl_gauss(uint64 stat_sum_1, uint64 stat_sum_x, uint64 stat_sum_x2,
 376                                                                                    uint64 n_dim, uint64 SUM_X, uint64 SUM_X2){
 377   uint64 stat_sum_xe ;  /* somme des xn region exterieure */
 378   uint32 ne ;             /* nombre de pixel region exterieure */
 379   double sigi2, sige2; /* variance region interieure et exterieure */
 380
 381   /* variance des valeurs des niveaux de gris a l'interieur du snake */
 382   sigi2 =
 383     ((double)stat_sum_x2/(double)stat_sum_1) -
 384     ((double)stat_sum_x/(uint64)stat_sum_1)*((double)stat_sum_x/(uint64)stat_sum_1) ;
 385
 386   /* variance des valeurs des niveaux de gris a l'exterieur du snake */
 387   ne = n_dim-stat_sum_1 ;
 388   stat_sum_xe = SUM_X - stat_sum_x ;
 389   sige2 =
 390     ((double)SUM_X2-stat_sum_x2)/(double)ne -
 391     ((double)stat_sum_xe/(uint64)ne)*((double)stat_sum_xe/(uint64)ne) ;
 392
 393   if ((sigi2 > 0)|(sige2 > 0))
 394   return  0.5*((double)stat_sum_1*log(sigi2) + (double)ne*log(sige2)) ;
 395   return -1 ;
 396 }
 397
 398
 399 __global__ void calcul_stats_snake(snake_node_gpu * d_snake, int  nnodes, int64 * d_stats_snake, double * vrais_min,
 400                                                                    t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2, int * TABLE_CODAGE, uint32 l
 401                                                                    )
 402 {
 403
 404   int id_nx, id_nprec, id_nprecprec ;
 405   int code_noeud, code_segment, pos ;
 406   __shared__ int64 s_stats_snake[3] ;
 407
 408   //init stats en shared mem
 409   s_stats_snake[0] = 0 ;
 410   s_stats_snake[1] = 0 ;
 411   s_stats_snake[2] = 0 ;
 412
 413
 414   for (id_nx = 0; id_nx < nnodes; id_nx++)
 415     {
 416           if (id_nx == 0) id_nprec = nnodes - 1;
 417           else id_nprec = id_nx - 1;
 418           if (id_nprec == 0) id_nprecprec = nnodes -1 ;
 419           else id_nprecprec = id_nprec - 1 ;
 420       /* gestion des segments partant du noeud */
 421       /* vers le noeud suivant dans l'ordre trigo */
 422       code_segment = d_snake[id_nprec].code_segment ;
 423       if (code_segment > 0)
 424                 {
 425                   /* on somme les contributions */
 426                   s_stats_snake[0] += d_snake[id_nprec].sum_1 ;
 427                   s_stats_snake[1] += d_snake[id_nprec].sum_x ;
 428                   s_stats_snake[2] += d_snake[id_nprec].sum_x2 ;
 429                 }
 430       else if (code_segment < 0)
 431                 {
 432                   /* on soustrait les contributions */
 433                   s_stats_snake[0] -= d_snake[id_nprec].sum_1 ;
 434                   s_stats_snake[1] -= d_snake[id_nprec].sum_x ;
 435                   s_stats_snake[2] -= d_snake[id_nprec].sum_x2 ;
 436                 }
 437       // else (code_segment == 0), on ne fait rien
 438       /* gestion des pixels connectant les segments */
 439       /* pixel de depart du segment actuel np --> np->noeud_suiv */
 440       /* freeman_out = np->freeman_out ; */
 441       /* freeman_in = np->noeud_prec->freeman_in ; */
 442           pos = d_snake[id_nprecprec].freeman_in*8 + d_snake[id_nprec].freeman_out ;
 443       code_noeud = TABLE_CODAGE[pos] ;
 444           pos = d_snake[id_nprec].posi*l + d_snake[id_nprec].posj ;
 445
 446       if (code_noeud > 0)
 447                 {
 448                   /* on somme les contributions */
 449                   s_stats_snake[0] += 1 + d_snake[id_nprec].posj ;
 450                   s_stats_snake[1] += cumul_x[pos] ;
 451                   s_stats_snake[2] += cumul_x2[pos] ;
 452                 }
 453       else if (code_noeud < 0)
 454                 {
 455                   /* on soustrait les contributions */
 456                   s_stats_snake[0] -= 1 + d_snake[id_nprec].posj ;
 457                   s_stats_snake[1] -= cumul_x[pos] ;
 458                   s_stats_snake[2] -= cumul_x2[pos] ;
 459                 }
 460       // else (code_pixel == 0), on ne fait rien
 461     }
 462   d_stats_snake[0] = s_stats_snake[0] ;
 463   d_stats_snake[1] = s_stats_snake[1] ;
 464   d_stats_snake[2] = s_stats_snake[2] ;
 465
 466   *vrais_min = codage_gl_gauss(s_stats_snake[0], s_stats_snake[1], s_stats_snake[2],
 467                                                            d_stats_snake[3], d_stats_snake[4], d_stats_snake[5]);
 468 }