src/lib_kernels_contribs.cu

   1 #include "constantes.h"
   2
   3 /*
   4  * determine et retourne le nb de pixels composant le segment (i,j)-(i1,j1)
   5  *
   6  */
   7 __device__ int calcul_nb_pixels(int i, int j, int i1, int j1){
   8            int absi, absj,nbpix =0;
   9            //MAX( ABS(i1-i) , ABS(j1-j)) + 1
  10            if (i1 > i) absi = i1 - i ; else absi = i - i1 ;
  11            if (j1 > j) absj = j1 - j ; else absj = j - j1 ;
  12            if (absi > absj ) nbpix = absi+1 ; else nbpix = absj+1 ;
  13            return nbpix ;
  14 }
  15
  16 /*
  17   construit la liste des coordonnées des 8  positions a tester pour l'ensemble des N noeuds du snake
  18   ainsi que le liste des nombres de pixels correspondant (en z, seg avant, en w seg apres)
  19   a executer avec N blocs de 8 threads
  20 */
  21 __global__ void liste_positions_a_tester(snake_node_gpu * d_snake, uint4 * liste_positions, uint32 * nb_pix_max, int pas, int nb_nodes, int h, int l){
  22   int tib = threadIdx.x;                                    // une position par thread
  23   int node = blockIdx.x;                                    // 1 noeud par bloc de threads
  24   int node_prec, node_suiv ;                                // indices des nodes
  25   int i, j, i1, j1, i2, j2, i3, j3, npix_prec, npix_suiv ;  // coordonnees des noeuds , nb de pixels
  26
  27   //lecture des coordonnees
  28   node_prec = (node == 0)? (nb_nodes-1) : (node - 1) ;
  29   node_suiv = (node == nb_nodes-1)? (0) : (node + 1) ;
  30   i1 = d_snake[node_prec].posi ;
  31   j1 = d_snake[node_prec].posj ;
  32   i2 = d_snake[node].posi ;
  33   j2 = d_snake[node].posj ;
  34   i3 = d_snake[node_suiv].posi ;
  35   j3 = d_snake[node_suiv].posj ;
  36
  37   switch(tib){                         // on considere un voisinage a 8 points
  38   case 0:
  39         i = i2 ;
  40         if ((j2 + pas) < l) j = j2 + pas ; else j = l-2 ;
  41         break;
  42   case 1:
  43         if ((j2 + pas) < l) j = j2 + pas ; else j = l-2 ;
  44         if (i2 > pas ) i = i2 - pas ; else i = 1 ;
  45         break;
  46   case 2:
  47         if (i2 > pas ) i = i2 - pas ; else i = 1 ;
  48         j = j2 ;
  49         break;
  50   case 3:
  51         if (i2 > pas) i = i2 - pas ; else i = 1 ;
  52         if (j2 > pas) j = j2 - pas ; else j = 1 ;
  53         break;
  54   case 4:
  55         i = i2 ;
  56         if (j2 > pas) j = j2 - pas ; else j = 1 ;
  57         break;
  58   case 5:
  59         if ((i2 + pas) < h) i = i2 + pas ; else i = h-2 ;
  60         if (j2 > pas) j = j2 - pas ; else j = 1 ;
  61         break;
  62   case 6:
  63         if ((i2 + pas) < h) i = i2 + pas ; else i = h-2 ;
  64         j = j2 ;
  65         break;
  66   case 7:
  67         if ((i2 + pas) < h) i = i2 + pas ; else i = h-2 ;
  68         if ((j2 + pas) < l) j = j2 + pas ; else j = l-2 ;
  69         break;
  70   }
  71
  72   //calcul des nombres de pixels dans chaque cas
  73   npix_prec =  calcul_nb_pixels(i,j,i1,j1) ;
  74   npix_suiv =  calcul_nb_pixels(i,j,i3,j3) ;
  75
  76   liste_positions[8*node + tib] = make_uint4(i,j, (uint)npix_prec, (uint)npix_suiv );
  77
  78   // calcule du maximum global du nb de pixels
  79   if ((node == 0) && (tib == 0))
  80         {
  81           *nb_pix_max = 0 ;
  82           for (int n=0; n<nb_nodes; n++)
  83                 {
  84                   if (liste_positions[n].z > *nb_pix_max) *nb_pix_max = liste_positions[n].z ;
  85                   if (liste_positions[n].w > *nb_pix_max) *nb_pix_max = liste_positions[n].w ;
  86                 }
  87         }
  88 }
  89
  90 /*
  91  * calcule :
  92  *  - les coordonnees de chaque pixel de chaque segment a evaluer pour les 8 voisins de chaque noeud (pair ou impair)
  93  *    cela represente 16 segments par noeud pair/impair
  94  *  - les contributions de chacun de ces pixels
  95  *  - les sommes, par blocs, des contributions des pixels (les sommes totales sont faites par le kernel somsom)
  96  *  - le code de chaque segment envisage
  97  */
  98
  99 __global__ void calcul_contribs_segments_blocs_full(snake_node_gpu * d_snake, int nb_nodes, uint4 * liste_points, uint32 npix_max,
 100                                                                                                         t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2,  int * d_codes_x16,
 101                                                                                                         int l, uint2 * liste_pix, uint64 * gsombloc, bool pairs)
 102 {
 103   // indices des elements
 104   int blockSize = blockDim.x ;                                                          // nb threads par bloc
 105   int tib = threadIdx.x ;                                                                       // position du thread dans le bloc
 106   int nblocs_noeud =  gridDim.x / nb_nodes ;                            // nb de blocs dédié à chaque noeud
 107   int nblocs_seg = nblocs_noeud / 16 ;                                  // nb de blocs dédiés à un segment de test
 108   int idx = blockDim.x*blockIdx.x + threadIdx.x ;                       // position absolue du thread ds la grille
 109   int id_interval = blockIdx.x / nblocs_noeud ;                         // indice de l'intervalle du noeud dans la grille
 110   int segment = ( blockIdx.x - id_interval*nblocs_noeud )/nblocs_seg ;  // indice du segment de 0 à 15
 111   int tis = idx - ( id_interval*nblocs_noeud + segment*nblocs_seg )*blockDim.x ; // position du thread ds le segment
 112   int id_base_seg = 16*id_interval + segment ;                          // indice du segment courant
 113   int id_base_pix = 5*id_base_seg ;                                                 // indice du pixel 0/5 du segment courant dans la liste_pix pour le calcul des freemans
 114
 115   //tab pour contribs pixels
 116   extern __shared__ tcontribs scumuls[];
 117
 118   //coordonnees des extremites de segment
 119   uint x1, y1, x2, y2 ;
 120   //indices des noeuds precedent(n1), courant(n2), suivant(n3)
 121   int n1, n2, n3 ;
 122   // pixel courant
 123   uint2 p ;
 124   // nb de pixels du segment precedent, suivant
 125   int xprec, xsuiv ;
 126
 127   // determine les indices des noeuds prec, courant, suiv
 128   n1 = id_interval -1 ;
 129   n2 = id_interval ;
 130   n3 = id_interval +1 ;
 131
 132   //gestion du bouclage du snake
 133   if (n1 < 0) n1 = nb_nodes-1 ;
 134   if (n3 >= nb_nodes) n3 = 0 ;
 135
 136
 137   //affectation des differentes positions aux différents segments 'blocs de threads'
 138   if ( segment < 8 ){
 139         x1 = d_snake[n1].posj ;
 140         y1 = d_snake[n1].posi ;
 141         x2 = liste_points[8*n2 + segment].y ;
 142         y2 = liste_points[8*n2 + segment].x ;
 143   } else {
 144         x1 = liste_points[8*n2 + segment-8].y ;
 145         y1 = liste_points[8*n2 + segment-8].x ;
 146         x2 = d_snake[n3].posj ;
 147         y2 = d_snake[n3].posi ;
 148   }
 149
 150   //params des deplacements
 151   int dx=x2-x1;
 152   int dy=y2-y1;
 153   uint abs_dx = ABS(dx);
 154   uint abs_dy = ABS(dy);
 155   uint nb_pix = abs_dy>abs_dx?(abs_dy+1):(abs_dx+1);
 156   int incx=0, incy=0;
 157
 158
 159   //calcul liste des pixels du segment (x1,y1)-(x2,y2)
 160   if (dy > 0) incy=1; else incy=-1 ;
 161   if (dx > 0) incx=1; else incx=-1 ;
 162
 163   if (tis < nb_pix){
 164         if (abs_dy > abs_dx){
 165           //1 thread par ligne pente 1/k
 166           double k = (double)dx/dy ;
 167           //coordonnees pixel
 168           p.x = y1 + incy*tis ;
 169           p.y = x1 + floor((double)incy*k*tis+0.5) ;
 170         } else {
 171           //1 thread par colonne pente k
 172           double k=(double)dy/dx ;
 173           //coordonnees pixel
 174           p.x = y1 + floor((double)(incx*k*tis)+0.5) ;
 175           p.y =  x1 + incx*tis ;
 176           // ordonnees des pixels suivant & precedent
 177           if ( tis > 0 ){
 178                 xsuiv = y1 + floor((double)(incx*k*(tis+1))+0.5) ;
 179                 xprec = y1 + floor((double)(incx*k*(tis-1))+0.5) ;
 180           }
 181         }
 182         // memorisation des valeurs necessaires au calcul des freemans et des centres
 183         if (tis == 0)  liste_pix[id_base_pix] = make_uint2(p.x, p.y) ;
 184         if (tis == 1)  liste_pix[id_base_pix +1] = make_uint2(p.x,p.y) ;
 185         if (tis == nb_pix/2) liste_pix[id_base_pix +2] = make_uint2(p.x,p.y) ;
 186         if (tis == nb_pix-2) liste_pix[id_base_pix +3] = make_uint2(p.x,p.y) ;
 187         if (tis == nb_pix-1) liste_pix[id_base_pix +4] = make_uint2(p.x,p.y) ;
 188   }
 189
 190   __syncthreads();
 191   // calcul contribs individuelles des pixels
 192   // dans le cas des segments 'plutot horizontaux', on ne garde qu'une contrib par ligne ; celle du pixel le plus a l'interieur du snake
 193   if ( (tis >0) && (tis < nb_pix-1)
 194            && ( ((abs_dy <= abs_dx) && ( ( xprec > p.x) || ( xsuiv > p.x)))
 195                         || (abs_dy > abs_dx) ) )
 196         {
 197           int pos = p.x * l + p.y ;
 198           scumuls[ CFI(tib)].c1 = 1 + p.y ;
 199           scumuls[ CFI(tib)].cx = cumul_x[ pos ] ;
 200           scumuls[CFI(tib)].cx2 = cumul_x2[ pos ];
 201         } else {
 202           scumuls[ CFI(tib)].c1 = 0;
 203           scumuls[ CFI(tib)].cx = 0;
 204           scumuls[ CFI(tib)].cx2 = 0;
 205   }
 206
 207    __syncthreads();
 208   // somme des contribs individuelles
 209   // unroll des  sommes partielles en shared memory
 210   if (blockSize >= 512) {
 211         if (tib < 256) {
 212           scumuls[ CFI(tib)].c1 += scumuls[ CFI(tib + 256) ].c1;
 213           scumuls[ CFI(tib)].cx += scumuls[ CFI(tib + 256) ].cx;
 214           scumuls[ CFI(tib)].cx2 += scumuls[CFI(tib + 256) ].cx2;
 215         }
 216         __syncthreads();
 217   }
 218
 219   if (blockSize >= 256) {
 220         if (tib < 128) {
 221           scumuls[ CFI(tib)].c1 += scumuls[ CFI(tib + 128) ].c1;
 222           scumuls[ CFI(tib)].cx += scumuls[ CFI(tib + 128) ].cx;
 223           scumuls[ CFI(tib)].cx2 += scumuls[CFI(tib + 128) ].cx2;
 224         }
 225         __syncthreads();
 226   }
 227   if (blockSize >= 128) {
 228         if (tib <  64) {
 229           scumuls[ CFI(tib)].c1 += scumuls[ CFI(tib +  64) ].c1;
 230           scumuls[ CFI(tib)].cx += scumuls[ CFI(tib +  64) ].cx;
 231           scumuls[CFI(tib)].cx2 += scumuls[CFI(tib +  64) ].cx2;
 232         }
 233         __syncthreads();
 234   }
 235
 236   //32 threads <==> 1 warp
 237  volatile tcontribs * scontribs = scumuls ;
 238  if (tib < 32){
 239          {
 240                 scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib + 32) ].c1;
 241                 scontribs[ CFI(tib)].cx += scontribs[ CFI(tib + 32) ].cx;
 242                 scontribs[CFI(tib)].cx2 += scontribs[CFI(tib + 32) ].cx2;
 243           }
 244  if (tib < 16)
 245           {
 246                 scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib + 16) ].c1;
 247                 scontribs[ CFI(tib)].cx += scontribs[ CFI(tib + 16) ].cx;
 248                 scontribs[CFI(tib)].cx2 += scontribs[CFI(tib + 16) ].cx2;
 249           }
 250  if (tib < 8)
 251           {
 252                 scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  8) ].c1;
 253                 scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  8) ].cx;
 254                 scontribs[CFI(tib)].cx2 += scontribs[CFI(tib +  8) ].cx2;
 255           }
 256  if (tib < 4)
 257         {
 258           scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  4) ].c1;
 259           scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  4) ].cx;
 260           scontribs[ CFI(tib)].cx2 +=scontribs[ CFI(tib +  4) ].cx2;
 261         }
 262  if (tib < 2)
 263         {
 264           scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  2) ].c1;
 265           scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  2) ].cx;
 266           scontribs[CFI(tib)].cx2 += scontribs[CFI(tib +  2) ].cx2;
 267          }
 268  if (tib == 0)
 269         {
 270           scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  1) ].c1;
 271           scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  1) ].cx;
 272           scontribs[CFI(tib)].cx2 += scontribs[CFI(tib +  1) ].cx2;
 273
 274   // resultat sommes partielles en gmem
 275   //if (tib == 0) {
 276         gsombloc[ blockIdx.x ] = scontribs[0].c1 ;
 277         gsombloc[ blockIdx.x + gridDim.x ] = scontribs[0].cx ;
 278         gsombloc[ blockIdx.x + 2*gridDim.x ] = scontribs[0].cx2 ;
 279
 280  //code seg
 281         if (dy > 0 ) d_codes_x16[id_base_seg] = -1 ;
 282         if (dy < 0 ) d_codes_x16[id_base_seg] =  1 ;
 283         if (dy == 0) d_codes_x16[id_base_seg]=  0 ;
 284   }
 285  }
 286 }
 287
 288
 289
 290 /*
 291  calcul des freeman et du centre de chaque segment de test
 292  a executer sur 'n_interval' blocs de 16 threads
 293  soit un thread par segment
 294 */
 295 __global__ void calcul_freemans_centre(uint2 * liste_pix, int * d_table_freeman, uint4 * d_freemans_x16){
 296
 297   int id_segment = threadIdx.x ;
 298   int id_freeman = ( blockDim.x*blockIdx.x + id_segment ) ;
 299   int id_base_pix = 5*id_freeman ;
 300
 301   //calculs freemans, centre et code segment
 302   //1 uint4 par segment
 303   int Dio, Djo, Dii, Dji;
 304
 305         //freeman out
 306         Dio = 1 + liste_pix[id_base_pix +1].x - liste_pix[id_base_pix].x ;
 307         Djo = 1 + liste_pix[id_base_pix +1].y - liste_pix[id_base_pix].y ;
 308         d_freemans_x16[id_freeman].z = d_table_freeman[3*Dio + Djo] ;
 309
 310
 311         //freeman_in
 312         Dii = 1 + liste_pix[id_base_pix +4].x - liste_pix[id_base_pix +3].x ;
 313         Dji = 1 + liste_pix[id_base_pix +4].y - liste_pix[id_base_pix +3].y ;
 314         d_freemans_x16[id_freeman].w = d_table_freeman[3*Dii + Dji] ;
 315
 316         //centre
 317         d_freemans_x16[id_freeman].x = liste_pix[id_base_pix +2].x ;
 318         d_freemans_x16[id_freeman].y = liste_pix[id_base_pix +2].y ;
 319
 320 }
 321
 322
 323 /*
 324   calcul des contribs 1, x et x2 des 16 segments
 325   autour du noeud
 326   a partir des sommes partielles des blocs
 327   1 bloc / 1 thread par segment
 328 */
 329
 330 __global__ void somsom_full(uint64 * somblocs, int nb_nodes, unsigned int nb_bl_seg, uint64 * somsom){
 331
 332   uint64 sdata[3];
 333   unsigned int seg = blockIdx.x ;
 334   unsigned int nb_seg = gridDim.x ;
 335
 336   //un thread par segment
 337   {
 338         sdata[0] = 0;
 339         sdata[1] = 0;
 340         sdata[2] = 0;
 341   }
 342
 343   for (int b=0; b < nb_bl_seg ; b++){ //voir atomicadd 64bits
 344         sdata[0] += somblocs[seg*nb_bl_seg + b];
 345         sdata[1] += somblocs[(seg + nb_seg)*nb_bl_seg + b];
 346         sdata[2] += somblocs[(seg + 2*nb_seg)*nb_bl_seg + b];
 347   }
 348
 349   //totaux en gmem
 350   {
 351         somsom[3*seg] = sdata[0];
 352         somsom[3*seg + 1] = sdata[1];
 353         somsom[3*seg + 2] = sdata[2];
 354   }
 355 }
 356
 357 /*
 358   version GPU de la fonction definie ds src/lib_math.c
 359 */
 360 __device__ bool test_inf_gpu(double arg1, double arg2){
 361    if (arg2 > 0)
 362          return (arg1 < (arg2*COEF_DECROI)) ;
 363   else
 364     return (arg1 < (arg2*INV_COEF_DECROI)) ;
 365 }
 366
 367 /*
 368   version GPU de la fonction codage_niveau_gris_hyp_gaussienne
 369  */
 370 __device__ double codage_gl_hyp_gauss(uint64 stat_sum_1, uint64 stat_sum_x, uint64 stat_sum_x2,
 371                                                                           uint64 n_dim, uint64 SUM_X, uint64 SUM_X2){
 372   uint64 stat_sum_xe ;  /* somme des xn region exterieure */
 373   uint32 ne ;             /* nombre de pixel region exterieure */
 374   double sigi2, sige2; /* variance region interieure et exterieure */
 375
 376   /* variance des valeurs des niveaux de gris a l'interieur du snake */
 377   sigi2 =
 378     (double)stat_sum_x2/(double)stat_sum_1 -
 379     ((double)stat_sum_x/stat_sum_1)*((double)stat_sum_x/stat_sum_1) ;
 380
 381   /* variance des valeurs des niveaux de gris a l'exterieur du snake */
 382   ne = n_dim-stat_sum_1 ;
 383   stat_sum_xe = SUM_X - stat_sum_x ;
 384   sige2 =
 385     ((double)SUM_X2-stat_sum_x2)/ne -
 386     ((double)stat_sum_xe/ne)*((double)stat_sum_xe/ne) ;
 387
 388   if ((sigi2 > 0)&&(sige2 > 0))
 389   return  0.5*((double)stat_sum_1*log(sigi2) + (double)ne*log(sige2)) ;
 390   return -1 ;
 391 }
 392
 393 /*
 394   soustrait, pour chaque intervalle [N1--Nx--N2]
 395   les contributions des segments N1--Nx et Nx--N2
 396
 397   A executer par nnodes blocs de 1 thread par intervalle
 398  */
 399
 400 __global__ void soustrait_aux_stats_2N_segments_noeud(snake_node_gpu * d_snake, int64 * d_stats_snake, int64 * d_stats_ref,
 401                                                                                                           t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2,
 402                                                                                                           int * TABLE_CODAGE, uint32 l
 403                                                                                                           )
 404 {
 405   int nnodes = gridDim.x ;
 406   int id_nx, id_nprec, id_nprecprec, id_nsuiv ;
 407   int code_noeud, pos;
 408   __shared__ int64 s_stats_snake[3] ;
 409
 410   id_nx = blockIdx.x ;
 411   if (id_nx == 0) id_nprec = nnodes - 1;
 412   else id_nprec = id_nx - 1;
 413   if (id_nprec == 0) id_nprecprec = nnodes -1 ;
 414   else id_nprecprec = id_nprec - 1 ;
 415   if (id_nx == nnodes-1) id_nsuiv = 0;
 416   else id_nsuiv = id_nx + 1 ;
 417
 418   //init avec les valeurs du contour actuel
 419   s_stats_snake[0] = d_stats_snake[0] ;
 420   s_stats_snake[1] = d_stats_snake[1] ;
 421   s_stats_snake[2] = d_stats_snake[2] ;
 422
 423   /* segment Na -- Nx */
 424   if (d_snake[id_nprec].code_segment > 0)
 425     {
 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 (d_snake[id_nprec].code_segment < 0)
 431     {
 432          s_stats_snake[0] += d_snake[id_nprec].sum_1 ;
 433      s_stats_snake[1] += d_snake[id_nprec].sum_x ;
 434      s_stats_snake[2] += d_snake[id_nprec].sum_x2 ;
 435     }
 436   // else (code_segment_NaNx == 0), on ne fait rien
 437
 438   /* gestion des pixels connectant les segments */
 439   /* pixel de depart du segment Na --> Nx */
 440   pos = d_snake[id_nprecprec].freeman_in*8 + d_snake[id_nprec].freeman_out ;
 441   code_noeud = TABLE_CODAGE[pos] ;
 442   pos = d_snake[id_nprec].posi*l + d_snake[id_nprec].posj ;
 443   if (code_noeud > 0)
 444     {
 445      s_stats_snake[0] -= 1 + d_snake[id_nprec].posj ;
 446      s_stats_snake[1] -= cumul_x[pos] ;
 447          s_stats_snake[2] -= cumul_x2[pos];
 448     }
 449   else if (code_noeud < 0)
 450     {
 451          s_stats_snake[0] += 1 + d_snake[id_nprec].posj ;
 452      s_stats_snake[1] += cumul_x[pos] ;
 453          s_stats_snake[2] += cumul_x2[pos];
 454     }
 455   // else (code_noeud == 0), on ne fait rien
 456
 457   /* -------------------------- */
 458   /* -------------------------- */
 459   /* segment Nx -- Nb */
 460   if ( d_snake[id_nx].code_segment > 0 )
 461     {
 462       // on soustrait la contribution
 463          s_stats_snake[0] -= d_snake[id_nx].sum_1 ;
 464      s_stats_snake[1] -= d_snake[id_nx].sum_x ;
 465      s_stats_snake[2] -= d_snake[id_nx].sum_x2 ;
 466     }
 467   else if ( d_snake[id_nx].code_segment < 0)
 468     {
 469          s_stats_snake[0] += d_snake[id_nx].sum_1 ;
 470      s_stats_snake[1] += d_snake[id_nx].sum_x ;
 471      s_stats_snake[2] += d_snake[id_nx].sum_x2 ;
 472     }
 473   // else (code_segment_NxNb == 0), on ne fait rien
 474
 475   /* gestion des pixels connectant les segments */
 476   /* pixel de depart du segment Nx --> Nb */
 477   pos = d_snake[id_nprec].freeman_in*8 + d_snake[id_nx].freeman_out ;
 478   code_noeud = TABLE_CODAGE[pos] ;
 479   pos = d_snake[id_nx].posi*l + d_snake[id_nx].posj ;
 480   if (code_noeud > 0)
 481     {
 482      s_stats_snake[0] -= 1 + d_snake[id_nx].posj  ;
 483      s_stats_snake[1] -= cumul_x[pos] ;
 484          s_stats_snake[2] -= cumul_x2[pos];
 485     }
 486   else if (code_noeud < 0)
 487     {
 488          s_stats_snake[0] += 1 + d_snake[id_nx].posj ;
 489      s_stats_snake[1] += cumul_x[pos] ;
 490          s_stats_snake[2] += cumul_x2[pos];
 491     }
 492   // else (code_noeud == 0), on ne fait rien
 493
 494   /* pixel d'arrivee du segment Nx --> Nb */
 495   pos = d_snake[id_nx].freeman_in*8 + d_snake[id_nsuiv].freeman_out ;
 496   code_noeud = TABLE_CODAGE[pos] ;
 497   pos = d_snake[id_nsuiv].posi*l +  d_snake[id_nsuiv].posj ;
 498   if (code_noeud > 0)
 499     {
 500          s_stats_snake[0] -= 1 + d_snake[id_nsuiv].posj ;
 501      s_stats_snake[1] -= cumul_x[pos] ;
 502          s_stats_snake[2] -= cumul_x2[pos];
 503     }
 504   else if (code_noeud < 0)
 505     {
 506          s_stats_snake[0] += 1 + d_snake[id_nsuiv].posj ;
 507      s_stats_snake[1] += cumul_x[pos] ;
 508          s_stats_snake[2] += cumul_x2[pos];
 509          }
 510   // else (code_noeud == 0), on ne fait rien
 511   __syncthreads();
 512
 513   d_stats_ref[3*id_nx] = s_stats_snake[0] ;
 514   d_stats_ref[3*id_nx + 1] = s_stats_snake[1] ;
 515   d_stats_ref[3*id_nx + 2] = s_stats_snake[2] ;
 516
 517 }
 518
 519
 520 /*
 521   calcul des stats associees a chaque position de test
 522
 523   EXEC : sur n_interval blocs de 8 threads
 524 */
 525 __global__ void calcul_stats_full(snake_node_gpu * d_snake, snake_node_gpu * d_snake_tmp, int nnodes, bool pairs, int64 * d_stats_snake,
 526                                                                   int64 * d_stats_ref, int64 * d_stats, uint64 * d_contribs,
 527                                                                   uint4 * d_liste_points, int * code_segment, uint4 * d_freemans,
 528                                                                   int * d_table_codes, t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2,
 529                                                                   uint32 h, uint32 l, double * vrais, double * vrais_min, bool * move){
 530
 531   int interval = blockIdx.x ;
 532   int seg = threadIdx.x ;
 533   int thread = seg;
 534   int id_nx, id_nprec, id_nprecprec, id_nsuiv, id_seg ;
 535   int code_noeud;
 536   __shared__ int64 s_stats_snake[3*8] ;
 537
 538   id_nx = interval ;
 539   if (id_nx == 0) id_nprec = nnodes - 1 ;
 540   else id_nprec = id_nx - 1 ;
 541   if (id_nprec == 0) id_nprecprec = nnodes -1 ;
 542   else id_nprecprec = id_nprec - 1 ;
 543   if (id_nx == nnodes-1) id_nsuiv = 0 ;
 544   else id_nsuiv = id_nx + 1 ;
 545
 546   //chargement en smem , prevoir CFI car on a 24x64bits par blocs => conflits
 547   s_stats_snake[3*thread + 0] = d_stats_ref[3*id_nx] ;
 548   s_stats_snake[3*thread + 1] = d_stats_ref[3*id_nx + 1] ;
 549   s_stats_snake[3*thread + 2] = d_stats_ref[3*id_nx + 2] ;
 550
 551
 552   //stats segments N1-Nx
 553   id_seg = 16*interval + seg ;
 554   if ( code_segment[id_seg] > 0 ){
 555         s_stats_snake[3*thread +0]  += d_contribs[3*id_seg] ;
 556         s_stats_snake[3*thread +1]  += d_contribs[3*id_seg + 1 ] ;
 557         s_stats_snake[3*thread +2]  += d_contribs[3*id_seg + 2 ] ;
 558   } else if ( code_segment[16*interval + seg] < 0 ) {
 559         s_stats_snake[3*thread +0]  -= d_contribs[3*id_seg] ;
 560         s_stats_snake[3*thread +1]  -= d_contribs[3*id_seg + 1] ;
 561         s_stats_snake[3*thread +2]  -= d_contribs[3*id_seg + 2] ;
 562   }
 563
 564   //stats noeud N1(i1,j1)
 565   int fo_N1 = d_freemans[id_seg].z ;
 566   int fi_Nprecprec = d_snake[id_nprecprec].freeman_in ;
 567   int pos = d_snake[id_nprec].posi*l + d_snake[id_nprec].posj ;
 568
 569   code_noeud = d_table_codes[fi_Nprecprec*8 + fo_N1];
 570   if (code_noeud > 0){
 571         s_stats_snake[3*thread +0] += 1 + d_snake[id_nprec].posj ;
 572         s_stats_snake[3*thread +1] += cumul_x[ pos ] ;
 573         s_stats_snake[3*thread +2] += cumul_x2[pos ] ;
 574   } else if (code_noeud < 0){
 575         s_stats_snake[3*thread +0] -= 1 + d_snake[id_nprec].posj ;
 576         s_stats_snake[3*thread +1] -= cumul_x[ pos ] ;
 577         s_stats_snake[3*thread +2] -= cumul_x2[pos ] ;
 578   }
 579
 580   //stats noeud Nx
 581   int fo_Nx = d_freemans[id_seg + 8].z ;
 582   int fi_Nx = d_freemans[id_seg].w ;
 583   int Nxi = d_liste_points[8*id_nx + seg].x ;
 584   int Nxj = d_liste_points[8*id_nx + seg].y ;
 585   pos = Nxi*l + Nxj ;
 586
 587   code_noeud = d_table_codes[fi_Nx*8 + fo_Nx];
 588   if (code_noeud > 0){
 589         s_stats_snake[3*thread +0] += 1 + Nxj ;
 590         s_stats_snake[3*thread +1] += cumul_x[ pos ] ;
 591         s_stats_snake[3*thread +2] += cumul_x2[pos ] ;
 592   }
 593   if (code_noeud < 0){
 594         s_stats_snake[3*thread +0] -= 1 + Nxj ;
 595         s_stats_snake[3*thread +1] -= cumul_x[ pos ] ;
 596         s_stats_snake[3*thread +2] -= cumul_x2[pos ] ;
 597   }
 598
 599   //stats segments Nx-N2
 600   seg += 8;
 601   id_seg = 16*interval + seg ;
 602   if ( code_segment[id_seg] > 0 ){
 603         s_stats_snake[3*thread +0]  += d_contribs[3*id_seg] ;
 604         s_stats_snake[3*thread +1]  += d_contribs[3*id_seg + 1] ;
 605         s_stats_snake[3*thread +2]  += d_contribs[3*id_seg + 2] ;
 606   }
 607   if ( code_segment[id_seg] < 0 ) {
 608         s_stats_snake[3*thread +0]  -= d_contribs[3*id_seg] ;
 609         s_stats_snake[3*thread +1]  -= d_contribs[3*id_seg + 1] ;
 610         s_stats_snake[3*thread +2]  -= d_contribs[3*id_seg + 2] ;
 611   }
 612
 613   //stats noeud N2(i2,j2)
 614   int fi_N2 = d_freemans[id_seg].w ;
 615   int fo_N2 = d_snake[id_nsuiv].freeman_out ;
 616   pos = d_snake[id_nsuiv].posi*l + d_snake[id_nsuiv].posj ;
 617
 618   code_noeud = d_table_codes[fi_N2*8 + fo_N2];
 619   if (code_noeud > 0){
 620         s_stats_snake[3*thread +0] += 1 + d_snake[id_nsuiv].posj ;
 621         s_stats_snake[3*thread +1] += cumul_x[ pos ] ;
 622         s_stats_snake[3*thread +2] += cumul_x2[pos ] ;
 623   }
 624   if (code_noeud < 0){
 625         s_stats_snake[3*thread +0] -= 1 + d_snake[id_nsuiv].posj ;
 626         s_stats_snake[3*thread +1] -= cumul_x[ pos ] ;
 627         s_stats_snake[3*thread +2] -= cumul_x2[pos ] ;
 628   }
 629
 630           //TODO
 631           //voir si on peut s'en passer
 632   d_stats[3*(8*interval + thread)] = s_stats_snake[3*thread +0];
 633   d_stats[3*(8*interval + thread) + 1] = s_stats_snake[3*thread +1];
 634   d_stats[3*(8*interval + thread) + 2] = s_stats_snake[3*thread +2];
 635
 636   //codage hyp gaussienne
 637   uint64 stat_sum_xe[8] ;    //somme des xn region exterieure
 638   uint32 ne[8] ;             // nombre de pixels region exterieure
 639   double sigi2[8], sige2[8]; //  carres des variances, regions interieure et exterieure
 640
 641   /* variance des valeurs des niveaux de gris a l'interieur du snake */
 642   sigi2[thread] =
 643         ((double)s_stats_snake[3*thread +2]/(double)s_stats_snake[3*thread +0]) -
 644         ((double)s_stats_snake[3*thread +1]/s_stats_snake[3*thread +0])*((double)s_stats_snake[3*thread +1]/s_stats_snake[3*thread +0]) ;
 645
 646   /* variance des valeurs des niveaux de gris a l'exterieur du snake */
 647   ne[thread] = h*l-s_stats_snake[3*thread +0] ;
 648   stat_sum_xe[thread] = d_stats_snake[4] - s_stats_snake[3*thread +1] ;
 649   sige2[thread] =
 650         (double)(d_stats_snake[5]-s_stats_snake[3*thread +2])/(double)ne[thread] -
 651         ((double)stat_sum_xe[thread]/ne[thread])*((double)stat_sum_xe[thread]/ne[thread]) ;
 652
 653   if (sige2[thread]>0 && sigi2[thread]>0)
 654         vrais[8*interval + thread] = 0.5*((double)s_stats_snake[3*thread]*log(sigi2[thread]) + (double)ne[thread]*log(sige2[thread])) ;
 655   else
 656         vrais[8*interval + thread] = -1.0;
 657
 658   if ( thread == 0 ){
 659          //init move
 660         move[id_nx] = false ;
 661         int pos_optim = -1;
 662         double vrais_tmp = *vrais_min;
 663         for (int v=0; v < 8; v++){
 664           if ( (vrais[8*interval + v] > 0) && (vrais[8*interval + v] < vrais_tmp*COEF_DECROI) ) {
 665                 vrais_tmp = vrais[8*interval + v];
 666                 pos_optim = v;
 667           }
 668         }
 669         if ( (pos_optim >-1) && !croisement(d_snake, id_nx, d_liste_points[8*id_nx + pos_optim].x, d_liste_points[8*id_nx + pos_optim].y, nnodes) )
 670                 {
 671                   /*maj data snake*/
 672                   move[id_nx] = true ;
 673                   //new position
 674                   d_snake_tmp[id_nx].posi = d_liste_points[8*id_nx + pos_optim].x ;
 675                   d_snake_tmp[id_nx].posj = d_liste_points[8*id_nx + pos_optim].y ;
 676                 }
 677         else {
 678                   //keep old position
 679                   d_snake_tmp[id_nx].posi = d_snake[id_nx].posi ;
 680                   d_snake_tmp[id_nx].posj = d_snake[id_nx].posj ;
 681
 682         }
 683   }
 684
 685
 686 }
 687
 688 __global__ void recalcul_stats_snake(snake_node_gpu * d_snake, int  nnodes, int64 * d_stats_snake, double * vrais_min,
 689                                                                          t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2, int * TABLE_CODAGE, uint32 l
 690                                                                          )
 691 {
 692
 693   int id_nx, id_nprec, id_nprecprec ;
 694   int code_noeud, code_segment, pos;
 695   int64 s_stats_snake[3] ;
 696
 697   //init stats en shared mem
 698   s_stats_snake[0] = 0 ;
 699   s_stats_snake[1] = 0 ;
 700   s_stats_snake[2] = 0 ;
 701
 702
 703   for (id_nx = 0; id_nx < nnodes; id_nx++)
 704     {
 705           if (id_nx == 0) id_nprec = nnodes - 1;
 706           else id_nprec = id_nx - 1;
 707           if (id_nprec == 0) id_nprecprec = nnodes -1 ;
 708           else id_nprecprec = id_nprec - 1 ;
 709       /* gestion des segments partant du noeud */
 710       /* vers le noeud suivant dans l'ordre trigo */
 711       code_segment = d_snake[id_nprec].code_segment ;
 712       if (code_segment > 0)
 713                 {
 714                   /* on somme les contributions */
 715                   s_stats_snake[0] += d_snake[id_nprec].sum_1 ;
 716                   s_stats_snake[1] += d_snake[id_nprec].sum_x ;
 717                   s_stats_snake[2] += d_snake[id_nprec].sum_x2 ;
 718                 }
 719       else if (code_segment < 0)
 720                 {
 721                   /* on soustrait les contributions */
 722                   s_stats_snake[0] -= d_snake[id_nprec].sum_1 ;
 723                   s_stats_snake[1] -= d_snake[id_nprec].sum_x ;
 724                   s_stats_snake[2] -= d_snake[id_nprec].sum_x2 ;
 725                 }
 726       // else (code_segment == 0), on ne fait rien
 727       /* gestion des pixels connectant les segments */
 728           pos = d_snake[id_nprecprec].freeman_in*8 + d_snake[id_nprec].freeman_out ;
 729       code_noeud = TABLE_CODAGE[pos] ;
 730           pos = d_snake[id_nprec].posi*l + d_snake[id_nprec].posj ;
 731
 732       if (code_noeud > 0)
 733                 {
 734                   /* on somme les contributions */
 735                   s_stats_snake[0] += 1 + d_snake[id_nprec].posj ;
 736                   s_stats_snake[1] += cumul_x[pos] ;
 737                   s_stats_snake[2] += cumul_x2[pos] ;
 738                 }
 739       else if (code_noeud < 0)
 740                 {
 741                   /* on soustrait les contributions */
 742                   s_stats_snake[0] -= 1 + d_snake[id_nprec].posj ;
 743                   s_stats_snake[1] -= cumul_x[pos] ;
 744                   s_stats_snake[2] -= cumul_x2[pos] ;
 745                 }
 746       // else (code_pixel == 0), on ne fait rien
 747     }
 748   d_stats_snake[0] = s_stats_snake[0] ;
 749   d_stats_snake[1] = s_stats_snake[1] ;
 750   d_stats_snake[2] = s_stats_snake[2] ;
 751
 752   *vrais_min = codage_gl_hyp_gauss(s_stats_snake[0], s_stats_snake[1], s_stats_snake[2],
 753                                                                    d_stats_snake[3], d_stats_snake[4], d_stats_snake[5]);
 754 }
 755
 756
 757 __global__ void ajoute_noeuds(snake_node_gpu * snake, snake_node_gpu * snake_tmp, int nnodes, int seuil, int * new_nb_nodes){
 758
 759
 760  int id_cpy = 0;
 761  for (int id_nx=0; id_nx < nnodes; id_nx++){
 762         //position du noeud existant
 763         snake_tmp[id_cpy].posi =  snake[id_nx].posi ;
 764     snake_tmp[id_cpy].posj =  snake[id_nx].posj ;
 765
 766         id_cpy++ ;
 767
 768         if ( snake_tmp[id_nx].nb_pixels > seuil)
 769         {
 770                 //position du nouveau noeud
 771                 snake_tmp[id_cpy].posi = snake[id_nx].centre_i ;
 772                 snake_tmp[id_cpy].posj = snake[id_nx].centre_j ;
 773                 id_cpy++ ;
 774         }
 775  }
 776
 777  for (int id_nx=0; id_nx < id_cpy; id_nx++){
 778         snake[id_cpy].posi =  snake_tmp[id_nx].posi ;
 779     snake[id_cpy].posj =  snake_tmp[id_nx].posj ;
 780  }
 781
 782 *new_nb_nodes = id_cpy-nnodes ;
 783 }
 784
 785 // pour copier le contenu d'un snake dans un autre
 786 // execution sur nnodes blocs de 1 threads
 787 __global__ void copie_snake(snake_node_gpu * d_source, snake_node_gpu * d_dest){
 788   int id = blockIdx.x ;
 789   d_dest[id].posi =  d_source[id].posi ;
 790   d_dest[id].posj =  d_source[id].posj ;
 791 }
 792
 793
 794
 795 __global__ void recalcul_contribs_segments_snake(snake_node_gpu * d_snake, int nb_nodes,
 796                                                                                                  t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2,
 797                                                                                                  int l, uint2 * liste_pix, uint64 * gsombloc )
 798 {
 799   // indices des elements
 800   int blockSize = blockDim.x ;
 801   int tib = threadIdx.x ;
 802   int nblocs_seg =  gridDim.x / nb_nodes ;
 803   int idx = blockDim.x*blockIdx.x + threadIdx.x ;
 804   int segment = blockIdx.x / nblocs_seg ;
 805   int tis = idx - segment*nblocs_seg*blockDim.x ;
 806
 807   //tab pour contribs pixels
 808   extern __shared__ tcontribs scumuls[];
 809
 810   //indices des noeuds
 811   uint x1, y1, x2, y2 ;
 812   int n1, n2 ;
 813   uint2 p ;
 814   int xsuiv, xprec ;
 815
 816   n1 = segment ;
 817   n2 = segment +1 ;
 818   //gestion du bouclage du snake
 819   if (n2 >= nb_nodes) n2 = 0 ;
 820
 821   //affectation des differentes positions aux différents segments 'blocs de threads'
 822         x1 = d_snake[n1].posj ;
 823         y1 = d_snake[n1].posi ;
 824         x2 = d_snake[n2].posj ;
 825         y2 = d_snake[n2].posi ;
 826
 827   //params des deplacements
 828   int dx=x2-x1;
 829   int dy=y2-y1;
 830   uint abs_dx = ABS(dx);
 831   uint abs_dy = ABS(dy);
 832   uint nb_pix = abs_dy>abs_dx?(abs_dy+1):(abs_dx+1); // alternative -> lecture ds liste_points[]
 833   int incx=0, incy=0;
 834
 835
 836   //calcul liste des pixels du segment (x1,y1)-(x2,y2)
 837   if (dy > 0) incy=1; else incy=-1 ;
 838   if (dx > 0) incx=1; else incx=-1 ;
 839
 840   if (tis < nb_pix){
 841         if (abs_dy > abs_dx){
 842           //1 thread par ligne
 843           double k = (double)dx/dy ;
 844           p.x = y1 + incy*tis ;
 845           p.y = x1 + floor((double)incy*k*tis+0.5) ;
 846
 847         } else {
 848           //1 thread par colonne
 849           double k=(double)dy/dx ;
 850           p.x = y1 + floor((double)(incx*k*tis)+0.5) ;
 851           p.y =  x1 + incx*tis ;
 852           if ( tis > 0 ){
 853                 xsuiv = y1 + floor((double)(incx*k*(tis+1))+0.5) ;
 854                 xprec = y1 + floor((double)(incx*k*(tis-1))+0.5) ;
 855           }
 856
 857         }
 858         if (tis == 0)  liste_pix[5*segment] = make_uint2(p.x, p.y) ;
 859         if (tis == 1)  liste_pix[5*segment +1] = make_uint2(p.x,p.y) ;
 860     if (tis == nb_pix/2) liste_pix[5*segment +2] = make_uint2(p.x,p.y) ;
 861     if (tis == nb_pix-2) liste_pix[5*segment +3] = make_uint2(p.x,p.y) ;
 862     if (tis == nb_pix-1) liste_pix[5*segment +4] = make_uint2(p.x,p.y) ;
 863
 864   }
 865   __syncthreads();
 866
 867   //calcul contribs individuelles des pixels
 868
 869   if ( (tis >0) && (tis < nb_pix-1)
 870            && ( ( (abs_dy <= abs_dx) && ( ( xprec > p.x) || ( xsuiv > p.x)) )
 871                         || (abs_dy > abs_dx) ) )
 872         {
 873         int pos = p.x * l + p.y ;
 874         scumuls[ CFI(tib)].c1 = 1+p.y;
 875         scumuls[ CFI(tib)].cx = cumul_x[ pos ] ;
 876         scumuls[CFI(tib)].cx2 = cumul_x2[ pos ];
 877   } else {
 878         scumuls[ CFI(tib)].c1 = 0;
 879         scumuls[ CFI(tib)].cx = 0;
 880         scumuls[CFI(tib)].cx2 = 0;
 881   }
 882
 883    __syncthreads();
 884   //somme des contribs individuelles
 885   // unroll des  sommes partielles en smem
 886
 887   if (blockSize >= 512) {
 888         if (tib < 256) {
 889           scumuls[ CFI(tib)].c1 += scumuls[ CFI(tib + 256) ].c1;
 890           scumuls[ CFI(tib)].cx += scumuls[ CFI(tib + 256) ].cx;
 891           scumuls[CFI(tib)].cx2 += scumuls[CFI(tib + 256) ].cx2;
 892         }
 893         __syncthreads();
 894   }
 895
 896   if (blockSize >= 256) {
 897         if (tib < 128) {
 898           scumuls[ CFI(tib)].c1 += scumuls[ CFI(tib + 128) ].c1;
 899           scumuls[ CFI(tib)].cx += scumuls[ CFI(tib + 128) ].cx;
 900           scumuls[CFI(tib)].cx2 += scumuls[CFI(tib + 128) ].cx2;
 901         }
 902         __syncthreads();
 903   }
 904   if (blockSize >= 128) {
 905         if (tib <  64) {
 906           scumuls[ CFI(tib)].c1 += scumuls[ CFI(tib +  64) ].c1;
 907           scumuls[ CFI(tib)].cx += scumuls[ CFI(tib +  64) ].cx;
 908           scumuls[CFI(tib)].cx2 += scumuls[ CFI(tib +  64) ].cx2;
 909         }
 910         __syncthreads();
 911   }
 912
 913   //32 threads <==> 1 warp
 914   volatile tcontribs * scontribs = scumuls ;
 915   if (tib < 32)
 916         {
 917           {
 918                 scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib + 32) ].c1;
 919                 scontribs[ CFI(tib)].cx += scontribs[ CFI(tib + 32) ].cx;
 920                 scontribs[ CFI(tib)].cx2 +=scontribs[ CFI(tib + 32) ].cx2;
 921           }
 922           if (tib<16)
 923           {
 924                 scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib + 16) ].c1;
 925                 scontribs[ CFI(tib)].cx += scontribs[ CFI(tib + 16) ].cx;
 926                 scontribs[ CFI(tib)].cx2 +=scontribs[ CFI(tib + 16) ].cx2;
 927           }
 928           if (tib<8)
 929           {
 930                 scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  8) ].c1;
 931                 scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  8) ].cx;
 932                 scontribs[ CFI(tib)].cx2 +=scontribs[ CFI(tib +  8) ].cx2;}
 933         if (tib<4){
 934           scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  4) ].c1;
 935           scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  4) ].cx;
 936           scontribs[ CFI(tib)].cx2 +=scontribs[ CFI(tib +  4) ].cx2;
 937           }
 938           if (tib<2){
 939           scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  2) ].c1;
 940           scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  2) ].cx;
 941           scontribs[ CFI(tib)].cx2 +=scontribs[ CFI(tib +  2) ].cx2;
 942           }
 943           if (tib==0){
 944           scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  1) ].c1;
 945           scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  1) ].cx;
 946           scontribs[ CFI(tib)].cx2 +=scontribs[ CFI(tib +  1) ].cx2;
 947   // resultat sommes partielles en gmem
 948   //if (tib == 0) {
 949         gsombloc[ blockIdx.x ] = scontribs[0].c1;
 950         gsombloc[ blockIdx.x + gridDim.x ] = scontribs[0].cx;
 951         gsombloc[ blockIdx.x + 2*gridDim.x ] = scontribs[0].cx2;
 952
 953   //calculs code segment
 954         //code seg
 955         if (dy > 0 ) d_snake[segment].code_segment = -1 ;
 956         if (dy < 0 ) d_snake[segment].code_segment = 1 ;
 957         if (dy == 0) d_snake[segment].code_segment = 0 ;
 958   }
 959         }
 960 }
 961
 962 __global__ void recalcul_freemans_centre(snake_node_gpu * snake, uint2 * liste_pix, int * d_table_freeman){
 963
 964   int id_segment = blockIdx.x ;
 965   int nb_segments = gridDim.x ;
 966   int id_base_pix = 5*id_segment ;
 967
 968   //calculs freemans, centre et code segment
 969   //1 uint4 par segment
 970   int Dio, Djo, Dii, Dji;
 971
 972         //freeman out
 973         Dio = 1 + liste_pix[id_base_pix +1].x - liste_pix[id_base_pix].x ;
 974         Djo = 1 + liste_pix[id_base_pix +1].y - liste_pix[id_base_pix].y ;
 975         snake[id_segment].freeman_out = d_table_freeman[3*Dio + Djo] ;
 976
 977         //freeman_in
 978         Dii = 1 + liste_pix[id_base_pix +4].x - liste_pix[id_base_pix +3].x ;
 979         Dji = 1 + liste_pix[id_base_pix +4].y - liste_pix[id_base_pix +3].y ;
 980         snake[id_segment].freeman_in = d_table_freeman[3*Dii + Dji] ;
 981
 982         //centre
 983         snake[id_segment].centre_i = liste_pix[id_base_pix +2].x ;
 984         snake[id_segment].centre_j = liste_pix[id_base_pix +2].y ;
 985
 986         //nb_pixels
 987         int nd = id_segment ;
 988         int nd_suiv = nd + 1  ;
 989         if (nd_suiv >= nb_segments ) nd_suiv = 0 ;
 990         snake[id_segment].nb_pixels = calcul_nb_pixels(snake[nd].posi, snake[nd].posj, snake[nd_suiv].posi, snake[nd_suiv].posj) ;
 991 }
 992
 993 /*
 994   sommme des contribs par bloc -> contribs segment, pour le snake
 995
 996   execution sur : 1bloc / 1 thread par segment
 997  */
 998
 999 __global__ void resomsom_snake(uint64 * somblocs, int nb_nodes, unsigned int nb_bl_seg, snake_node_gpu * d_snake){
1000
1001   uint64 sdata[3];
1002   unsigned int seg = blockIdx.x ;
1003
1004   //un thread par segment
1005   {
1006         sdata[0] = 0;
1007         sdata[1] = 0;
1008         sdata[2] = 0;
1009   }
1010
1011   for (int b=0; b < nb_bl_seg ; b++){
1012         sdata[0] += somblocs[seg*nb_bl_seg + b];
1013         sdata[1] += somblocs[(seg + nb_nodes)*nb_bl_seg + b];
1014         sdata[2] += somblocs[(seg + 2*nb_nodes)*nb_bl_seg + b];
1015   }
1016
1017   //totaux en gmem
1018   {
1019         d_snake[seg].sum_1 = sdata[0];
1020         d_snake[seg].sum_x = sdata[1];
1021         d_snake[seg].sum_x2 = sdata[2];
1022   }
1023 }
1024