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 = 0;
  30                 d_snake[i].code_segment = 0;
  31
  32           }
  33   }
  34 }
  35
  36 __global__ void genere_snake_bande_gpu(snake_node_gpu * d_snake, int j1, int j2, int h){
  37   if (threadIdx.x == 0){
  38         int n = 0 ;
  39         /* n0 */
  40         d_snake[n].posi = 0 ;
  41         d_snake[n].posj = j1 ;
  42         n++ ;
  43         /* n1 */
  44         d_snake[n].posi = h-1 ;
  45         d_snake[n].posj = j1 ;
  46         n++ ;
  47         /* n2 */
  48         d_snake[n].posi = h-1 ;
  49         d_snake[n].posj = j2 ;
  50         n++ ;
  51         /* n3 */
  52         d_snake[n].posi = 0 ;
  53         d_snake[n].posj = j2 ;
  54
  55         for (int i=0; i<4; i++)
  56           {
  57                 d_snake[i].freeman_in = 0;
  58                 d_snake[i].freeman_out = 0;
  59                 d_snake[i].centre_i = 0;
  60                 d_snake[i].centre_j = 0;
  61                 d_snake[i].last_move = 0;
  62                 d_snake[i].nb_pixels = 0;
  63                 d_snake[i].code_segment = 0;
  64
  65           }
  66   }
  67 }
  68
  69
  70 __global__ void genere_snake_rectangle(snake_node_gpu * d_snake, int i1, int i2, int j1, int j2){
  71   if (threadIdx.x == 0){
  72         int n = 0 ;
  73         /* n0 */
  74         d_snake[n].posi = i1 ;
  75         d_snake[n].posj = j1 ;
  76         n++ ;
  77         /* n1 */
  78         d_snake[n].posi = i2 ;
  79         d_snake[n].posj = j1 ;
  80         n++ ;
  81         /* n2 */
  82         d_snake[n].posi = i2 ;
  83         d_snake[n].posj = j2 ;
  84         n++ ;
  85         /* n3 */
  86         d_snake[n].posi = i1 ;
  87         d_snake[n].posj = j2 ;
  88
  89         for (int i=0; i<4; i++)
  90           {
  91                 d_snake[i].freeman_in = 0;
  92                 d_snake[i].freeman_out = 0;
  93                 d_snake[i].centre_i = 0;
  94                 d_snake[i].centre_j = 0;
  95                 d_snake[i].last_move = 0;
  96                 d_snake[i].nb_pixels = 0;
  97                 d_snake[i].code_segment = 0;
  98
  99           }
 100   }
 101 }
 102
 103
 104 /*
 105   Evalue le critere dans l'hypothese gaussienne
 106  */
 107 __device__ double critere_gauss( int n, int ni, uint64 sxi, uint64 sxi2, uint64 sigX, uint64 sigX2){
 108   int ne = n - ni ;
 109   uint64 sxe = sigX - sxi ;
 110   uint64 sxe2= sigX2 - sxi2;
 111   double sigi2, sige2, critere ;
 112
 113   /* variance des valeurs des niveaux de gris a l'interieur */
 114   sigi2 =
 115         ((double)sxi2/ni) -
 116         ((double)sxi/ni)*((double)sxi/ni) ;
 117
 118   /* variance des valeurs des niveaux de gris a l'exterieur */
 119   sige2 =
 120         ((double)sxe2)/ne -
 121         ((double)sxe/ne)*((double)sxe/ne) ;
 122
 123   if ( (sigi2 > 0) && (sige2 > 0) )
 124         critere =  0.5*((double)ni*log(sigi2) + (double)ne*log(sige2)) ;
 125   else critere = DBL_MAX ;
 126   return critere;
 127 }
 128
 129
 130 __global__ void calcul_contribs_segments_snake(snake_node_gpu * d_snake, int nb_nodes,
 131                                                                                            t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2,
 132                                                                                            int l, uint2 * liste_pix, t_sum_x2 * gsombloc, int * d_table_freeman)
 133 {
 134   // indices des elements
 135   int blockSize = blockDim.x ;
 136   int tib = threadIdx.x ;
 137   int nblocs_seg =  gridDim.x / nb_nodes ;
 138   int idx = blockDim.x*blockIdx.x + threadIdx.x ;
 139   int segment = blockIdx.x / nblocs_seg ;
 140   int tis = idx - segment*nblocs_seg*blockDim.x ;
 141
 142   //tab pour coordonnées pixels & contribs pixels de taille = (blockDim.x+offset(dec,dec2) )*(sizeof(t_sum_1+t_sum_x+t_sum_x2))
 143   extern __shared__ t_sum_1 scumuls_1[] ; // blockDim varie selon la longueur des segments => taille smem dynamique
 144   t_sum_x* scumuls_x = (t_sum_x*) &scumuls_1[CFI(blockDim.x)] ;
 145   t_sum_x2* scumuls_x2 = (t_sum_x2*) &scumuls_x[CFI(blockDim.x)] ;
 146
 147   //indices des noeuds
 148   uint x1, y1, x2, y2 ;
 149   int n1, n2 ;
 150
 151   n1 = segment ;
 152   n2 = segment +1 ;
 153   //gestion du bouclage du snake
 154   if (n2 >= nb_nodes) n2 = 0 ;
 155
 156   //affectation des differentes positions aux différents segments 'blocs de threads'
 157         x1 = d_snake[n1].posj ;
 158         y1 = d_snake[n1].posi ;
 159         x2 = d_snake[n2].posj ;
 160         y2 = d_snake[n2].posi ;
 161
 162   //params des deplacements
 163   int dx=x2-x1;
 164   int dy=y2-y1;
 165   uint abs_dx = ABS(dx);
 166   uint abs_dy = ABS(dy);
 167   uint nb_pix = abs_dy>abs_dx?(abs_dy+1):(abs_dx+1); // alternative -> lecture ds liste_points[]
 168   int incx=0, incy=0;
 169   uint2 p ;
 170   int xprec, xsuiv ;
 171
 172   //calcul liste des pixels du segment (x1,y1)-(x2,y2)
 173   if (dy > 0) incy=1; else incy=-1 ;
 174   if (dx > 0) incx=1; else incx=-1 ;
 175
 176   if (tis < nb_pix){
 177         if (abs_dy > abs_dx){
 178           //1 thread par ligne
 179           double k = (double)dx/dy ;
 180           p.x = y1 + incy*tis ;
 181           p.y =  x1 + floor((double)incy*k*tis+0.5) ;
 182           //enreg. coords. pixels en global mem pour freemans
 183
 184           if ((tis < 2)||(tis > nb_pix - 3)||(tis == nb_pix/2))
 185                 {
 186                   liste_pix[idx].x = p.x ;
 187                   liste_pix[idx].y = p.y ;
 188                   }
 189
 190         } else {
 191           //1 thread par colonne
 192           double k=(double)dy/dx ;
 193           p.x = y1 + floor((double)(incx*k*tis)+0.5) ;
 194           p.y = x1 + incx*tis ;
 195           if ( tis > 0 ){
 196                 xsuiv = y1 + floor((double)(incx*k*(tis+1))+0.5) ;
 197                 xprec = y1 + floor((double)(incx*k*(tis-1))+0.5) ;
 198           }
 199           //enreg. coords. pixels en global mem pour freeman
 200           //TODO
 201           //on peut calculer les freemans des segments
 202           //sans stocker l'ensemble des valeurs des pixels
 203           //juste avec les derivees aux extremites a calculer ici
 204
 205           if ((tis < 2)||(tis > nb_pix - 3)||(tis == nb_pix/2))
 206                 {
 207                   liste_pix[idx].x = p.x ;
 208                   liste_pix[idx].y = p.y ;
 209                   }
 210
 211         }
 212   }
 213   __syncthreads();
 214
 215   //calcul contribs individuelles des pixels
 216
 217   if ( (tis >0) && (tis < nb_pix-1)
 218            && ( (abs_dy <= abs_dx)
 219                         && ( (xprec > p.x) || (xsuiv > p.x))
 220                         || (abs_dy > abs_dx) ) )
 221         {
 222         int pos = p.x * l + p.y ;
 223         scumuls_1[ CFI(tib)] = 1+p.y;
 224         scumuls_x[ CFI(tib)] = cumul_x[ pos ] ;
 225         scumuls_x2[CFI(tib)] = cumul_x2[ pos ];
 226   } else {
 227         scumuls_1[ CFI(tib)] = 0;
 228         scumuls_x[ CFI(tib)] = 0;
 229         scumuls_x2[CFI(tib)] = 0;
 230   }
 231
 232    __syncthreads();
 233   //somme des contribs individuelles
 234   // unroll des  sommes partielles en smem
 235
 236   if (blockSize >= 512) {
 237         if (tib < 256) {
 238           scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib + 256) ];
 239           scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib + 256) ];
 240           scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib + 256) ];
 241         }
 242         __syncthreads();
 243   }
 244
 245   if (blockSize >= 256) {
 246         if (tib < 128) {
 247           scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib + 128) ];
 248           scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib + 128) ];
 249           scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib + 128) ];
 250         }
 251         __syncthreads();
 252   }
 253   if (blockSize >= 128) {
 254         if (tib <  64) {
 255           scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib +  64) ];
 256           scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib +  64) ];
 257           scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib +  64) ];
 258         }
 259         __syncthreads();
 260   }
 261
 262   //32 threads <==> 1 warp
 263   if (tib < 32)
 264         {
 265           {
 266                 scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib + 32) ];
 267                 scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib + 32) ];
 268                 scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib + 32) ];
 269           }
 270           {
 271                 scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib + 16) ];
 272                 scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib + 16) ];
 273                 scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib + 16) ];
 274           }
 275           {
 276                 scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib +  8) ];
 277                 scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib +  8) ];
 278                 scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib +  8) ];
 279           }
 280           scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib +  4) ];
 281           scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib +  4) ];
 282           scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib +  4) ];
 283
 284           scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib +  2) ];
 285           scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib +  2) ];
 286           scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib +  2) ];
 287
 288           scumuls_1[ CFI(tib)] += scumuls_1[ CFI(tib +  1) ];
 289           scumuls_x[ CFI(tib)] += scumuls_x[ CFI(tib +  1) ];
 290           scumuls_x2[CFI(tib)] += scumuls_x2[CFI(tib +  1) ];
 291         }
 292
 293   // resultat sommes partielles en gmem
 294   if (tib == 0) {
 295         gsombloc[ blockIdx.x ] = (t_sum_x2) scumuls_1[0];
 296         gsombloc[ blockIdx.x + gridDim.x ] = (t_sum_x2) scumuls_x[0];
 297         gsombloc[ blockIdx.x + 2*gridDim.x ] = (t_sum_x2) scumuls_x2[0];
 298   }
 299
 300   //calculs freemans, centre et code segment
 301   //1 uint4 par segment
 302
 303   int Di, Dj;
 304   if (tis == 0){
 305         Di = 1 + liste_pix[idx+1].x - liste_pix[idx].x ;
 306         Dj = 1 + liste_pix[idx+1].y - liste_pix[idx].y ;
 307         d_snake[segment].freeman_out = d_table_freeman[3*Di + Dj] ;
 308         //code seg
 309         if (dy > 0 ) d_snake[segment].code_segment = -1 ;
 310         if (dy < 0 ) d_snake[segment].code_segment = 1 ;
 311         if (dy == 0) d_snake[segment].code_segment = 0 ;
 312   }
 313
 314   if (tis == nb_pix-1){
 315         Di = 1 + liste_pix[idx].x - liste_pix[idx-1].x  ;
 316         Dj = 1 + liste_pix[idx].y - liste_pix[idx-1].y;
 317         d_snake[segment].freeman_in = d_table_freeman[3*Di + Dj] ;
 318   }
 319
 320   if (tis == (nb_pix/2)){
 321         d_snake[segment].centre_i = liste_pix[idx].x ;
 322         d_snake[segment].centre_j = liste_pix[idx].y ;
 323         }
 324 }
 325
 326 /*
 327   sommme des contribs par bloc -> contribs segment, pour le snake
 328
 329   execution sur : 1bloc / 1 thread par segment
 330  */
 331
 332 __global__ void somsom_snake(t_sum_x2 * somblocs, int nb_nodes, unsigned int nb_bl_seg, snake_node_gpu * d_snake){
 333
 334   t_sum_x2 sdata[3];
 335   unsigned int seg = blockIdx.x ;
 336
 337   //un thread par segment
 338   {
 339         sdata[0] = 0;
 340         sdata[1] = 0;
 341         sdata[2] = 0;
 342   }
 343
 344   for (int b=0; b < nb_bl_seg ; b++){
 345         sdata[0] += somblocs[seg*nb_bl_seg + b];
 346         sdata[1] += somblocs[(seg + nb_nodes)*nb_bl_seg + b];
 347         sdata[2] += somblocs[(seg + 2*nb_nodes)*nb_bl_seg + b];
 348   }
 349
 350   //totaux en gmem
 351   {
 352         d_snake[seg].sum_1 = sdata[0];
 353         d_snake[seg].sum_x = sdata[1];
 354         d_snake[seg].sum_x2 = sdata[2];
 355   }
 356 }
 357
 358 __device__ double codage_gl_gauss(uint64 stat_sum_1, uint64 stat_sum_x, uint64 stat_sum_x2,
 359                                                                                    uint64 n_dim, uint64 SUM_X, uint64 SUM_X2){
 360   uint64 stat_sum_xe ;  /* somme des xn region exterieure */
 361   uint32 ne ;             /* nombre de pixel region exterieure */
 362   double sigi2, sige2; /* variance region interieure et exterieure */
 363
 364   /* variance des valeurs des niveaux de gris a l'interieur du snake */
 365   sigi2 =
 366     ((double)stat_sum_x2/(double)stat_sum_1) -
 367     ((double)stat_sum_x/(uint64)stat_sum_1)*((double)stat_sum_x/(uint64)stat_sum_1) ;
 368
 369   /* variance des valeurs des niveaux de gris a l'exterieur du snake */
 370   ne = n_dim-stat_sum_1 ;
 371   stat_sum_xe = SUM_X - stat_sum_x ;
 372   sige2 =
 373     ((double)SUM_X2-stat_sum_x2)/(double)ne -
 374     ((double)stat_sum_xe/(uint64)ne)*((double)stat_sum_xe/(uint64)ne) ;
 375
 376   if ((sigi2 > 0) && (sige2 > 0))
 377   return  0.5*((double)stat_sum_1*log(sigi2) + (double)ne*log(sige2)) ;
 378   return -1 ;
 379 }
 380
 381
 382 __global__ void calcul_stats_snake(snake_node_gpu * d_snake, int  nnodes, int64 * d_stats_snake, double * vrais_min,
 383                                                                    t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2, int * TABLE_CODAGE, uint32 l
 384                                                                    )
 385 {
 386
 387   int id_nx, id_nprec, id_nprecprec ;
 388   int code_noeud, code_segment, pos ;
 389   __shared__ int64 s_stats_snake[3] ;
 390
 391   //init stats en shared mem
 392   s_stats_snake[0] = 0 ;
 393   s_stats_snake[1] = 0 ;
 394   s_stats_snake[2] = 0 ;
 395
 396
 397   for (id_nx = 0; id_nx < nnodes; id_nx++)
 398     {
 399           if (id_nx == 0) id_nprec = nnodes - 1;
 400           else id_nprec = id_nx - 1;
 401           if (id_nprec == 0) id_nprecprec = nnodes -1 ;
 402           else id_nprecprec = id_nprec - 1 ;
 403       /* gestion des segments partant du noeud */
 404       /* vers le noeud suivant dans l'ordre trigo */
 405       code_segment = d_snake[id_nprec].code_segment ;
 406       if (code_segment > 0)
 407                 {
 408                   /* on somme les contributions */
 409                   s_stats_snake[0] += d_snake[id_nprec].sum_1 ;
 410                   s_stats_snake[1] += d_snake[id_nprec].sum_x ;
 411                   s_stats_snake[2] += d_snake[id_nprec].sum_x2 ;
 412                 }
 413       else if (code_segment < 0)
 414                 {
 415                   /* on soustrait les contributions */
 416                   s_stats_snake[0] -= d_snake[id_nprec].sum_1 ;
 417                   s_stats_snake[1] -= d_snake[id_nprec].sum_x ;
 418                   s_stats_snake[2] -= d_snake[id_nprec].sum_x2 ;
 419                 }
 420       // else (code_segment == 0), on ne fait rien
 421       /* gestion des pixels connectant les segments */
 422       /* pixel de depart du segment actuel np --> np->noeud_suiv */
 423       /* freeman_out = np->freeman_out ; */
 424       /* freeman_in = np->noeud_prec->freeman_in ; */
 425           pos = d_snake[id_nprecprec].freeman_in*8 + d_snake[id_nprec].freeman_out ;
 426       code_noeud = TABLE_CODAGE[pos] ;
 427           pos = d_snake[id_nprec].posi*l + d_snake[id_nprec].posj ;
 428
 429       if (code_noeud > 0)
 430                 {
 431                   /* on somme les contributions */
 432                   s_stats_snake[0] += 1 + d_snake[id_nprec].posj ;
 433                   s_stats_snake[1] += cumul_x[pos] ;
 434                   s_stats_snake[2] += cumul_x2[pos] ;
 435                 }
 436       else if (code_noeud < 0)
 437                 {
 438                   /* on soustrait les contributions */
 439                   s_stats_snake[0] -= 1 + d_snake[id_nprec].posj ;
 440                   s_stats_snake[1] -= cumul_x[pos] ;
 441                   s_stats_snake[2] -= cumul_x2[pos] ;
 442                 }
 443       // else (code_pixel == 0), on ne fait rien
 444     }
 445   d_stats_snake[0] = s_stats_snake[0] ;
 446   d_stats_snake[1] = s_stats_snake[1] ;
 447   d_stats_snake[2] = s_stats_snake[2] ;
 448
 449   *vrais_min = codage_gl_gauss(s_stats_snake[0], s_stats_snake[1], s_stats_snake[2],
 450                                                            d_stats_snake[3], d_stats_snake[4], d_stats_snake[5]);
 451 }
 452
 453 // kernel d'init rectangulaire au niveau snake
 454 // un block par point de base et un point opposé par thread du bloc
 455
 456 __global__ void calcul_contribs_snake4(snake_node_gpu * snake, t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2, int h, int l,
 457                                                                            int64 * sums, t_rectangle_snake * critere)
 458 {
 459   // nb de diagonales testees par bloc (ie. par point de base NO)
 460   int blockSize = blockDim.x ;
 461   // indice du second point de chaque diagonale (=Opposite Point, = point SE)
 462   int OPib = threadIdx.x ;
 463   // coordonnees de chaque point de base (NO)
 464   int BPi = blockIdx.x ;
 465   int BPj = blockIdx.y ;
 466   //coordonnees de chaque Opposite Point (SE)
 467   double incThread = ((h-BPi)*(l-BPj) + blockSize -1)/(double)blockSize ;
 468   int OPi = (double)(OPib*incThread)/(l - BPj) ;
 469   int OPj = OPib*incThread - OPi*(l-BPj) ;
 470   OPi += BPi ;
 471   OPj += BPj ;
 472   if (OPi >= h) OPi = h-1 ;
 473   if (OPj >= l) OPj = l-1 ;
 474   //indices des pixels dans les images cumulees
 475   int posG, posD;
 476   //contrib 1 du snake
 477   int C1 = (OPi - BPi)*(OPj - BPj) ;
 478
 479   //pour stocker contribs de chaque snake d'un block
 480   //TODO on peut utiliser une structure restreinte (sans le c1)
 481   extern __shared__ tcontribs scumuls[];
 482
 483   //calcul contribs du snake
 484   scumuls[CFI(OPib)].cx = 0;
 485   scumuls[CFI(OPib)].cx2 = 0;
 486   for (int k=BPi ; k < OPi ; k++)
 487         {
 488           posG = (BPi+k)*l + BPj ;
 489           posD = posG - BPj + OPj ;
 490           scumuls[CFI(OPib)].cx  += (cumul_x[ posD ] - cumul_x[ posG ]) ;
 491           scumuls[CFI(OPib)].cx2 += (cumul_x2[ posD ] - cumul_x2[ posG ]);
 492   }
 493
 494   //calcul de critère pour chaque snake
 495   uint64 stat_sum_xe ;  /* somme des xn region exterieure */
 496   uint32 ne ;           /* nombre de pixel region exterieure */
 497   double sigi2, sige2;  /* variance region interieure et exterieure */
 498   int index_crit;
 499
 500   /* variance des valeurs des niveaux de gris a l'interieur du snake */
 501   sigi2 =
 502     ((double)scumuls[CFI(OPib)].cx2/(double)C1) -
 503     ((double)scumuls[CFI(OPib)].cx/(uint64)C1)*((double)scumuls[CFI(OPib)].cx/(uint64)C1) ;
 504
 505   /* variance des valeurs des niveaux de gris a l'exterieur du snake */
 506
 507   ne = sums[3] - C1 ;
 508   stat_sum_xe = sums[4] - scumuls[CFI(OPib)].cx ;
 509   sige2 =
 510     ((double)sums[5]-scumuls[CFI(OPib)].cx2)/(double)ne -
 511     ((double)stat_sum_xe/(uint64)ne)*((double)stat_sum_xe/(uint64)ne) ;
 512   index_crit = blockSize*(BPi*gridDim.y + BPj) + OPib ;
 513
 514    critere[ index_crit ].bpi = BPi;
 515    critere[ index_crit ].bpj = BPj;
 516    critere[ index_crit ].opi = OPi;
 517    critere[ index_crit ].opj = OPj;
 518
 519   if ((sigi2 > 0)|(sige2 > 0))
 520         critere[ index_crit ].crit = 0.5*((double)C1*log(sigi2) + (double)ne*log(sige2)) ;
 521   else
 522         critere[ index_crit ].crit = -1 ;
 523
 524   // identification meilleur snake du bloc
 525   // laissé au CPU pour test mais le principe de ce kernel n'est pas efficace.
 526
 527 }
 528 __global__ void calcul_contribs_cols(t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2, int h, int l, tcontribs * sums )
 529 {
 530
 531   int bs  = blockDim.x ;       // nb pixels par bloc
 532   int nbC = gridDim.x ;        // nb de colonnes traitees
 533   int bpc = gridDim.y ;        // nb de blocs par colonne
 534   int idC = blockIdx.x ;       // indice de la colonne ds la grille
 535   int idB = blockIdx.y ;       // indice du bloc de pixels ds la colonne
 536   int tib = threadIdx.x ;      // indice du thread ds le bloc courant
 537
 538   int i = idB*bs + tib ;       // ligne du pixel
 539   int j = idC * (l/nbC) ;      // colonne du pixel
 540   int pos =  i*l + j ;         // index absolu du pixel ds l'image
 541
 542   extern __shared__ tcontribs scontrib[]; // pour stocker les contribs partielles
 543
 544   if (i < h)
 545         {
 546           scontrib[ CFI(tib) ].cx = cumul_x[ pos ] ;
 547           scontrib[ CFI(tib) ].cx2= cumul_x2[pos ] ;
 548         }
 549   else
 550         {
 551           scontrib[ CFI(tib) ].cx = 0 ;
 552           scontrib[ CFI(tib) ].cx2= 0 ;
 553         }
 554   __syncthreads();
 555
 556   // somme des contribs individuelles par bloc
 557   // unroll des  sommes partielles en smem
 558   if (bs >= 1024) {
 559         if (tib < 512) {
 560           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib + 512) ].cx;
 561           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib + 512) ].cx2;
 562         }
 563         __syncthreads();
 564   }
 565
 566   if (bs >= 512) {
 567         if (tib < 256) {
 568           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib + 256) ].cx;
 569           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib + 256) ].cx2;
 570         }
 571         __syncthreads();
 572   }
 573
 574   if (bs >= 256) {
 575         if (tib < 128) {
 576           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib + 128) ].cx;
 577           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib + 128) ].cx2;
 578         }
 579         __syncthreads();
 580   }
 581   if (bs >= 128) {
 582         if (tib <  64) {
 583           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib +  64) ].cx;
 584           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib +  64) ].cx2;
 585         }
 586         __syncthreads();
 587   }
 588
 589   //32 threads <==> 1 warp
 590   if (tib < 32)
 591         {
 592           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib + 32) ].cx;
 593           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib + 32) ].cx2;
 594         }
 595   if (tib < 16)
 596         {
 597           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib + 16) ].cx;
 598           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib + 16) ].cx2;
 599         }
 600   if (tib < 8)
 601         {
 602           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib +  8) ].cx;
 603           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib +  8) ].cx2;
 604         }
 605   if (tib < 4)
 606         {
 607           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib +  4) ].cx;
 608           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib +  4) ].cx2;
 609         }
 610   if (tib < 2)
 611         {
 612           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib +  2) ].cx;
 613           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib +  2) ].cx2;
 614         }
 615   if (tib < 1)
 616         {
 617           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib +  1) ].cx;
 618           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib +  1) ].cx2;
 619         }
 620
 621   __syncthreads();
 622
 623   // resultat contribs partielles ( par bloc) en gmem
 624   if (tib == 0) {
 625         // c1 s'obtient directement car les sommes sont à j constant ( j=idC )
 626         sums[ idC*bpc +idB ].cx = scontrib[0].cx ;
 627         sums[ idC*bpc +idB ].cx2  = scontrib[0].cx2 ;
 628   }
 629 }
 630
 631 /*
 632   sommme des contribs par bloc -> contribs colonnes
 633   execution sur : (1 bloc de 1 thread) par colonne
 634  */
 635
 636 __global__ void somsom_contribs(tcontribs * somblocs, int bpc, tcontribs * somcols){
 637
 638   tcontribs scontrib ;
 639   int col = blockIdx.x ;
 640   int base = col*bpc ;
 641
 642   //un thread par colonne
 643   {
 644         scontrib.cx = 0;
 645         scontrib.cx2 = 0;
 646   }
 647
 648   for (int b=0; b < bpc ; b++){
 649         scontrib.cx += somblocs[ base +b ].cx ;
 650         scontrib.cx2 += somblocs[ base +b ].cx2 ;
 651   }
 652
 653   //totaux en gmem
 654   {
 655         somcols[ col ].cx  = scontrib.cx ;
 656         somcols[ col ].cx2 = scontrib.cx2;
 657   }
 658
 659 }
 660
 661 /*
 662   recherche des criteres minis dans les blocs
 663  */
 664 __device__ void mini_critere_bloc(double2 *critere, int bs){
 665
 666   int tib = threadIdx.x ;
 667
 668   critere[ CFI(tib) ].y = tib ; //initialisation des indices pour les minima
 669
 670   if ( (bs >= 1024) && (tib < 512) ) {
 671         if ( critere[ CFI(tib + 512) ].x < critere[ CFI(tib) ].x ){
 672           critere[ CFI(tib) ] = critere[ CFI(tib + 512) ] ;
 673           __syncthreads() ;
 674         }
 675   }
 676
 677   if ( (bs >= 512) && (tib < 256) ) {
 678         if ( critere[ CFI(tib + 256) ].x < critere[ CFI(tib) ].x ){
 679           critere[ CFI(tib) ] = critere[ CFI(tib + 256) ] ;
 680           __syncthreads() ;
 681         }
 682   }
 683
 684   if ( (bs >= 256) && (tib < 128) ) {
 685         if ( critere[ CFI(tib + 128) ].x < critere[ CFI(tib) ].x ){
 686           critere[ CFI(tib) ] = critere[ CFI(tib + 128) ] ;
 687           __syncthreads() ;
 688         }
 689   }
 690
 691   if ( (bs >= 128) && (tib < 64) ) {
 692         if ( critere[ CFI(tib + 64) ].x < critere[ CFI(tib) ].x ){
 693           critere[ CFI(tib) ] = critere[ CFI(tib + 64) ] ;
 694           __syncthreads() ;
 695         }
 696   }
 697
 698   //32 threads <==> 1 warp
 699   if (tib < 32){
 700         if ( critere[ CFI(tib + 32) ].x < critere[ CFI(tib) ].x ){
 701           critere[ CFI(tib) ] = critere[ CFI(tib + 32) ] ;
 702         }
 703   }
 704
 705   if (tib < 16){
 706         if ( critere[ CFI(tib + 16) ].x < critere[ CFI(tib) ].x ){
 707           critere[ CFI(tib) ] = critere[ CFI(tib + 16) ] ;
 708         }
 709   }
 710
 711   if (tib < 8){
 712         if ( critere[ CFI(tib + 8) ].x < critere[ CFI(tib) ].x ){
 713           critere[ CFI(tib) ] = critere[ CFI(tib + 8) ] ;
 714         }
 715   }
 716
 717   if (tib < 4){
 718         if ( critere[ CFI(tib + 4) ].x < critere[ CFI(tib) ].x ){
 719           critere[ CFI(tib) ] = critere[ CFI(tib + 4) ] ;
 720         }
 721   }
 722
 723   if (tib < 2){
 724         if ( critere[ CFI(tib + 2) ].x < critere[ CFI(tib) ].x ){
 725           critere[ CFI(tib) ] = critere[ CFI(tib + 2) ] ;
 726         }
 727   }
 728
 729   if (tib < 1){
 730         if ( critere[ CFI(tib + 1) ].x < critere[ CFI(tib) ].x ){
 731           critere[ CFI(tib) ] = critere[ CFI(tib + 1) ] ;
 732         }
 733   }
 734   __syncthreads() ;
 735 }
 736
 737
 738 /*
 739   calcule, pour l'ensemble des permutations possibles
 740   des colonnes j1 et j2 (j1<J2) :
 741   - le critère correspondant
 742  */
 743 __global__ void calcul_contribs_permutations(tcontribs * somcols, int bpc, int h, int l, uint64 * stats_img, double2 * miniblocs){
 744
 745   int bs = blockDim.x ;
 746   int j1 = blockIdx.x ;
 747   int j2 = threadIdx.x ;
 748
 749   //shared mem dynamique 'critere' et 'indice'
 750   extern __shared__ double2 critere[] ;
 751
 752   critere[ CFI(j2) ].y = j2 ; //initialisation des indices pour les minima
 753
 754   __syncthreads() ;
 755
 756   if (j1 < j2){
 757         int colj1 = j1*bpc ;
 758         int colj2 = j2*bpc ;
 759         int ni = h*(colj2 - colj1) ;
 760         int ne = h*l - ni ;                                  // h*l ou stats_img[3]
 761         uint64 stat_sum_xi = somcols[j2].cx - somcols[j1].cx ;
 762         uint64 stat_sum_xi2 = somcols[j2].cx2 - somcols[j1].cx2 ;
 763         uint64 stat_sum_xe = stats_img[4] - stat_sum_xi ;    /* somme des xn region exterieure */
 764         uint64 stat_sum_xe2 = stats_img[5] - stat_sum_xi2 ;
 765     double sigi2, sige2;                                 /* variances regions interieure et exterieure */
 766
 767         /* variance des valeurs des niveaux de gris a l'interieur */
 768         sigi2 =
 769           ((double)stat_sum_xi2/ni) -
 770           ((double)stat_sum_xi/ni)*((double)stat_sum_xi/ni) ;
 771
 772         /* variance des valeurs des niveaux de gris a l'exterieur */
 773         sige2 =
 774           ((double)stat_sum_xe2)/ne -
 775           ((double)stat_sum_xe/ne)*((double)stat_sum_xe/ne) ;
 776
 777         if ((sigi2 > 0)|(sige2 > 0))
 778           critere[CFI(j2)].x =  0.5*((double)ni*log(sigi2) + (double)ne*log(sige2)) ;
 779         else critere[CFI(j2)].x = DBL_MAX;
 780   }
 781   else critere[CFI(j2)].x = DBL_MAX ;
 782
 783   __syncthreads(); // on s'assure que toutes les valeurs des criteres sont calculees avant de chercher le mini
 784
 785   /*
 786   if ( (bs >= 1024) && (j2 < 512) ) {
 787         if ( critere[ CFI(j2 + 512) ].x < critere[ CFI(j2) ].x ){
 788           critere[ CFI(j2) ] = critere[ CFI(j2 + 512) ] ;
 789           __syncthreads() ;
 790         }
 791   }
 792
 793   if ( (bs >= 512) && (j2 < 256) ) {
 794         if ( critere[ CFI(j2 + 256) ].x < critere[ CFI(j2) ].x ){
 795           critere[ CFI(j2) ] = critere[ CFI(j2 + 256) ] ;
 796           __syncthreads() ;
 797         }
 798   }
 799
 800   if ( (bs >= 256) && (j2 < 128) ) {
 801         if ( critere[ CFI(j2 + 128) ].x < critere[ CFI(j2) ].x ){
 802           critere[ CFI(j2) ] = critere[ CFI(j2 + 128) ] ;
 803           __syncthreads() ;
 804         }
 805   }
 806
 807   if ( (bs >= 128) && (j2 < 64) ) {
 808         if ( critere[ CFI(j2 + 64) ].x < critere[ CFI(j2) ].x ){
 809           critere[ CFI(j2) ] = critere[ CFI(j2 + 64) ] ;
 810           __syncthreads() ;
 811         }
 812   }
 813
 814   //32 threads <==> 1 warp
 815   if (j2 < 32){
 816         if ( critere[ CFI(j2 + 32) ].x < critere[ CFI(j2) ].x ){
 817           critere[ CFI(j2) ] = critere[ CFI(j2 + 32) ] ;
 818         }
 819   }
 820
 821   if (j2 < 16){
 822         if ( critere[ CFI(j2 + 16) ].x < critere[ CFI(j2) ].x ){
 823           critere[ CFI(j2) ] = critere[ CFI(j2 + 16) ] ;
 824         }
 825   }
 826
 827   if (j2 < 8){
 828         if ( critere[ CFI(j2 + 8) ].x < critere[ CFI(j2) ].x ){
 829           critere[ CFI(j2) ] = critere[ CFI(j2 + 8) ] ;
 830         }
 831   }
 832
 833   if (j2 < 4){
 834         if ( critere[ CFI(j2 + 4) ].x < critere[ CFI(j2) ].x ){
 835           critere[ CFI(j2) ] = critere[ CFI(j2 + 4) ] ;
 836         }
 837   }
 838
 839   if (j2 < 2){
 840         if ( critere[ CFI(j2 + 2) ].x < critere[ CFI(j2) ].x ){
 841           critere[ CFI(j2) ] = critere[ CFI(j2 + 2) ] ;
 842         }
 843   }
 844
 845   if (j2 < 1){
 846         if ( critere[ CFI(j2 + 1) ].x < critere[ CFI(j2) ].x ){
 847           critere[ CFI(j2) ] = critere[ CFI(j2 + 1) ] ;
 848         }
 849         //critere[0] contient maintenant les params du mini du bloc
 850         miniblocs[ j1 ] = critere[0] ;
 851         }
 852   */
 853   mini_critere_bloc(critere, bs) ;
 854   if (j2 == 0) miniblocs[ j1 ] = critere[0] ;
 855 }
 856
 857
 858 __device__ void somme_contribs_blocs(tcontribs * scontrib, int bs){
 859   int tib = threadIdx.x ;
 860   // somme des contribs individuelles par bloc
 861   // unroll des  sommes partielles en smem
 862   if (bs >= 1024) {
 863         if (tib < 512) {
 864           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib + 512) ].cx;
 865           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib + 512) ].cx2;
 866         }
 867         __syncthreads();
 868   }
 869
 870   if (bs >= 512) {
 871         if (tib < 256) {
 872           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib + 256) ].cx;
 873           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib + 256) ].cx2;
 874         }
 875         __syncthreads();
 876   }
 877
 878   if (bs >= 256) {
 879         if (tib < 128) {
 880           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib + 128) ].cx;
 881           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib + 128) ].cx2;
 882         }
 883         __syncthreads();
 884   }
 885   if (bs >= 128) {
 886         if (tib <  64) {
 887           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib +  64) ].cx;
 888           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib +  64) ].cx2;
 889         }
 890         __syncthreads();
 891   }
 892
 893   //32 threads <==> 1 warp
 894   if (tib < 32)
 895         {
 896           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib + 32) ].cx;
 897           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib + 32) ].cx2;
 898         }
 899   if (tib < 16)
 900         {
 901           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib + 16) ].cx;
 902           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib + 16) ].cx2;
 903         }
 904   if (tib < 8)
 905         {
 906           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib +  8) ].cx;
 907           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib +  8) ].cx2;
 908         }
 909   if (tib < 4)
 910         {
 911           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib +  4) ].cx;
 912           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib +  4) ].cx2;
 913         }
 914   if (tib < 2)
 915         {
 916           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib +  2) ].cx;
 917           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib +  2) ].cx2;
 918         }
 919   if (tib < 1)
 920         {
 921           scontrib[ CFI(tib) ].cx += scontrib[ CFI(tib +  1) ].cx;
 922           scontrib[ CFI(tib) ].cx2 += scontrib[CFI(tib +  1) ].cx2;
 923         }
 924
 925   __syncthreads();
 926 }
 927
 928 /*
 929   effectue la somme (pix par pix) ds un bloc des contribs sur les colonnes j1 et j2 (j1 < j2)
 930   Grille de calcul GPU :
 931   -> threads = bs selon capacités/perfs
 932   -> grid = div = (h+bs-1)/bs (h = nb lignes)
 933  */
 934 __global__ void calcul_contrib_conjuguee_colonnes(t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2, int h, int l, int j1, int j2, tcontribs * vect_contrib)
 935 {
 936   int bs = blockDim.x ;
 937   int id_bloc = blockIdx.x ;
 938   int tib = threadIdx.x ;
 939   int idx = id_bloc * bs + tib ;
 940
 941   extern __shared__ tcontribs sh_contrib[] ;
 942
 943   if (idx < h){
 944         sh_contrib[CFI(tib)].cx = cumul_x[ idx * l + j2 ] - cumul_x[ idx * l + j1 ] ;
 945         sh_contrib[CFI(tib)].cx2= cumul_x2[ idx * l + j2 ] - cumul_x2[ idx * l + j1 ] ;
 946   } else {
 947         sh_contrib[CFI(tib)].cx =  0 ;
 948         sh_contrib[CFI(tib)].cx2=  0 ;
 949   }
 950   __syncthreads(); // synchro avant de sommer
 951
 952   /*sommes par bloc*/
 953   somme_contribs_blocs(sh_contrib, bs);
 954
 955   /*enregistrement somme des blocs*/
 956   if (tib == 0) vect_contrib[ id_bloc ] = sh_contrib[0] ;
 957
 958 }
 959
 960 /*
 961   evalue les hauteurs les plus probables des horizontales
 962   du rectangle 'englobant' (qui ne l'est pas nécessairement en fait)
 963   Grille de calcul GPU :
 964   -> threads = div = nb de blocs de lignes (Cf. calcul_contrib_conjuguee_colonnes() )
 965   -> grid = div x div x 1
 966 */
 967 __global__ void calcul_contribs_snake_rectangle(tcontribs * vect_contribs, tcontribs * soms)
 968 {
 969   int bs = blockDim.x ;
 970   int div = gridDim.x ;
 971   int tib = threadIdx.x ;
 972   int i1 = blockIdx.x ;
 973   int i2 = blockIdx.y ;
 974
 975   extern __shared__ tcontribs scontribs[] ;
 976
 977   if (i2 >= i1){
 978         if ( (tib >= i1) && (tib <= i2) ) {
 979           scontribs[CFI(tib)].cx = vect_contribs[tib].cx ;
 980           scontribs[CFI(tib)].cx2= vect_contribs[tib].cx2;
 981         } else {
 982           scontribs[CFI(tib)].cx  = 0;
 983           scontribs[CFI(tib)].cx2 = 0;
 984         }
 985         __syncthreads() ;
 986         /*sommes par bloc*/
 987         somme_contribs_blocs(scontribs, bs);
 988
 989         /*ranger les sommes ds un vecteur en gmem*/
 990         if (tib == 0) soms[ i1*div + i2 ] = scontribs[0] ;
 991   } else {
 992         if (tib == 0) {
 993           soms[ i1*div + i2 ].cx = 0 ;
 994           soms[ i1*div + i2 ].cx = 0 ;
 995         }
 996   }
 997 }
 998
 999 /*
1000   Evalue la valeur du critere pour chaque combinaison possible
1001   des blocs de lignes définis par  calcul_contrib_conjuguee_colonnes()
1002   grille de calcul GPU :
1003   threads -> nextPow2 (div)
1004   grid -> div x 1 x1
1005
1006   TODO voir a inverser les indices i1/i2
1007  */
1008
1009 __global__ void calcul_critere_permutations_verticales(tcontribs *soms, int lpb, int j1, int j2, int h, int l, uint64 sigX, uint64 sigX2, double2 *miniblocs){
1010
1011   int n = h*l ;
1012   int bs = blockDim.x ;
1013   int div = gridDim.x ;
1014   int i2 = threadIdx.x ;
1015   int i1 = blockIdx.x ;
1016
1017   extern __shared__ double2 sh_mini[];
1018
1019   if ( (i2 < div) && (i2 >= i1) )
1020         sh_mini[ CFI(i2) ].x = critere_gauss(n, (i2-i1+1)*lpb*(j2 - j1), soms[ i1*div +i2 ].cx, soms[ i1*div + i2 ].cx2, sigX, sigX2);
1021   else sh_mini[ CFI(i2) ].x = DBL_MAX ;
1022   __syncthreads() ;
1023
1024   mini_critere_bloc(sh_mini, bs) ;
1025
1026   if (i2==0) {
1027         //miniblocs[i1] = sh_mini[0] ;
1028         miniblocs[i1].x = sh_mini[0].x ;
1029         miniblocs[i1].y = sh_mini[0].y ;
1030         }
1031
1032 }