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