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