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