#include "constantes.h"

/*
 * determine et retourne le nb de pixels composant le segment (i,j)-(i1,j1)
 *
 */
__device__ int calcul_nb_pixels(int i, int j, int i1, int j1){
	   int absi, absj,nbpix =0;
	   //MAX( ABS(i1-i) , ABS(j1-j)) + 1
	   if (i1 > i) absi = i1 - i ; else absi = i - i1 ;
	   if (j1 > j) absj = j1 - j ; else absj = j - j1 ;
	   if (absi > absj ) nbpix = absi+1 ; else nbpix = absj+1 ;
	   return nbpix ;
}

/*
  construit la liste des coordonnées des 8  positions a tester pour l'ensemble des N noeuds du snake
  ainsi que le liste des nombres de pixels correspondant (en z, seg avant, en w seg apres)
  a executer avec N blocs de 8 threads
*/
__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){
  int tib = threadIdx.x;                                    // une position par thread
  int node = blockIdx.x;                                    // 1 noeud par bloc de threads
  int node_prec, node_suiv ;                                // indices des nodes
  int i, j, i1, j1, i2, j2, i3, j3, npix_prec, npix_suiv ;  // coordonnees des noeuds , nb de pixels

  //lecture des coordonnees
  node_prec = (node == 0)? (nb_nodes-1) : (node - 1) ;
  node_suiv = (node == nb_nodes-1)? (0) : (node + 1) ;
  i1 = d_snake[node_prec].posi ;
  j1 = d_snake[node_prec].posj ;
  i2 = d_snake[node].posi ;
  j2 = d_snake[node].posj ;
  i3 = d_snake[node_suiv].posi ;
  j3 = d_snake[node_suiv].posj ;
  
  switch(tib){                         // on considere un voisinage a 8 points
  case 0:
	i = i2 ;
	if ((j2 + pas) < l) j = j2 + pas ; else j = l-2 ; 
	break;
  case 1:
	if ((j2 + pas) < l) j = j2 + pas ; else j = l-2 ;
	if (i2 > pas ) i = i2 - pas ; else i = 1 ;
	break;
  case 2:
	if (i2 > pas ) i = i2 - pas ; else i = 1 ;
	j = j2 ; 
	break;
  case 3:
	if (i2 > pas) i = i2 - pas ; else i = 1 ;
	if (j2 > pas) j = j2 - pas ; else j = 1 ;
	break;
  case 4:
	i = i2 ;
	if (j2 > pas) j = j2 - pas ; else j = 1 ;
	break;
  case 5:
	if ((i2 + pas) < h) i = i2 + pas ; else i = h-2 ;
	if (j2 > pas) j = j2 - pas ; else j = 1 ;
	break;
  case 6:
	if ((i2 + pas) < h) i = i2 + pas ; else i = h-2 ;
	j = j2 ;
	break;
  case 7:
	if ((i2 + pas) < h) i = i2 + pas ; else i = h-2 ;
	if ((j2 + pas) < l) j = j2 + pas ; else j = l-2 ;
	break;
  }

  //calcul des nombres de pixels dans chaque cas
  npix_prec =  calcul_nb_pixels(i,j,i1,j1) ;
  npix_suiv =  calcul_nb_pixels(i,j,i3,j3) ;
  
  liste_positions[8*node + tib] = make_uint4(i,j, (uint)npix_prec, (uint)npix_suiv );
  
  // calcule du maximum global du nb de pixels
  if ((node == 0) && (tib == 0))
	{
	  *nb_pix_max = 0 ;
	  for (int n=0; n<nb_nodes; n++)
		{
		  if (liste_positions[n].z > *nb_pix_max) *nb_pix_max = liste_positions[n].z ;
		  if (liste_positions[n].w > *nb_pix_max) *nb_pix_max = liste_positions[n].w ;
		}
	}
}

/*
 * calcule :
 *  - les coordonnees de chaque pixel de chaque segment a evaluer pour les 8 voisins de chaque noeud (pair ou impair)
 *    cela represente 16 segments par noeud pair/impair
 *  - les contributions de chacun de ces pixels
 *  - les sommes, par blocs, des contributions des pixels (les sommes totales sont faites par le kernel somsom)
 *  - le code de chaque segment envisage
 */

__global__ void calcul_contribs_segments_blocs_full(snake_node_gpu * d_snake, int nb_nodes, uint4 * liste_points, uint32 npix_max,
													t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2,  int * d_codes_x16,
													int l, uint2 * liste_pix, uint64 * gsombloc, bool pairs)
{
  // indices des elements
  int blockSize = blockDim.x ;					                        // nb threads par bloc
  int tib = threadIdx.x ;						                        // position du thread dans le bloc
  int nblocs_noeud =  gridDim.x / (nb_nodes/2 + pairs*(nb_nodes%2)) ;   // nb de blocs dédié à chaque noeud 
  int nblocs_seg = nblocs_noeud / 16 ;                                  // nb de blocs dédiés à un segment de test
  int idx = blockDim.x*blockIdx.x + threadIdx.x ;                       // position absolue du thread ds la grille
  int id_interval = blockIdx.x / nblocs_noeud ;                         // indice de l'intervalle du noeud dans la grille
  int segment = ( blockIdx.x - id_interval*nblocs_noeud )/nblocs_seg ;  // indice du segment de 0 à 15
  int tis = idx - ( id_interval*nblocs_noeud + segment*nblocs_seg )*blockDim.x ; // position du thread ds le segment
  int id_base_seg = 16*id_interval + segment ;                          // indice du segment courant
  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

  //tab pour contribs pixels 
  extern __shared__ tcontribs scumuls[];  

  //coordonnees des extremites de segment
  uint x1, y1, x2, y2 ;
  //indices des noeuds precedent(n1), courant(n2), suivant(n3)
  int n1, n2, n3 ;
  // pixel courant
  uint2 p ;
  // nb de pixels du segment precedent, suivant
  int xprec, xsuiv ; 

  // determine les indices des noeuds prec, courant, suiv
  if (pairs)
	{
	  n1 = 2*id_interval -1 ;
	  n2 = 2*id_interval ;
	  n3 = 2*id_interval +1 ;
	}
  else
	{
	  n1 = 2*id_interval ;
	  n2 = 2*id_interval +1 ;
	  n3 = 2*id_interval +2 ;
	}
  //gestion du bouclage du snake
  if (n1 < 0) n1 = nb_nodes-1 ;
  if (n3 >= nb_nodes) n3 = 0 ;

  
  //affectation des differentes positions aux différents segments 'blocs de threads'
  if ( segment < 8 ){
	x1 = d_snake[n1].posj ;
	y1 = d_snake[n1].posi ; 
	x2 = liste_points[8*n2 + segment].y ;
	y2 = liste_points[8*n2 + segment].x ;
  } else {
	x1 = liste_points[8*n2 + segment-8].y ;
	y1 = liste_points[8*n2 + segment-8].x ;
	x2 = d_snake[n3].posj ;
	y2 = d_snake[n3].posi ;
  }
  
  //params des deplacements
  int dx=x2-x1;
  int dy=y2-y1;
  uint abs_dx = ABS(dx);
  uint abs_dy = ABS(dy);
  uint nb_pix = abs_dy>abs_dx?(abs_dy+1):(abs_dx+1);
  int incx=0, incy=0;
  
  
  //calcul liste des pixels du segment (x1,y1)-(x2,y2)  
  if (dy > 0) incy=1; else incy=-1 ;
  if (dx > 0) incx=1; else incx=-1 ;
  
  if (tis < nb_pix){
	if (abs_dy > abs_dx){
	  //1 thread par ligne pente 1/k
	  double k = (double)dx/dy ;
	  //coordonnees pixel
	  p.x = y1 + incy*tis ;
	  p.y = x1 + floor((double)incy*k*tis+0.5) ;
	} else {
	  //1 thread par colonne pente k	  
	  double k=(double)dy/dx ;
	  //coordonnees pixel
	  p.x = y1 + floor((double)(incx*k*tis)+0.5) ;
	  p.y =  x1 + incx*tis ;
	  // ordonnees des pixels suivant & precedent
	  if ( tis > 0 ){ 
		xsuiv = y1 + floor((double)(incx*k*(tis+1))+0.5) ;
		xprec = y1 + floor((double)(incx*k*(tis-1))+0.5) ;
	  }
	}
	// memorisation des valeurs necessaires au calcul des freemans et des centres
	if (tis == 0)  liste_pix[id_base_pix] = make_uint2(p.x, p.y) ;
	if (tis == 1)  liste_pix[id_base_pix +1] = make_uint2(p.x,p.y) ;
	if (tis == nb_pix/2) liste_pix[id_base_pix +2] = make_uint2(p.x,p.y) ;
	if (tis == nb_pix-2) liste_pix[id_base_pix +3] = make_uint2(p.x,p.y) ;
	if (tis == nb_pix-1) liste_pix[id_base_pix +4] = make_uint2(p.x,p.y) ;
  }
  
  __syncthreads();  
  // calcul contribs individuelles des pixels
  // 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
  if ( (tis >0) && (tis < nb_pix-1)
	   && ( ((abs_dy <= abs_dx) && ( ( xprec > p.x) || ( xsuiv > p.x)))
			|| (abs_dy > abs_dx) ) )
	{
	  int pos = p.x * l + p.y ;
	  scumuls[ CFI(tib)].c1 = 1 + p.y ;  
	  scumuls[ CFI(tib)].cx = cumul_x[ pos ] ;
	  scumuls[CFI(tib)].cx2 = cumul_x2[ pos ];
	} else {
	  scumuls[ CFI(tib)].c1 = 0;
	  scumuls[ CFI(tib)].cx = 0;
	  scumuls[ CFI(tib)].cx2 = 0;
  }
  
   __syncthreads();
  // somme des contribs individuelles 
  // unroll des  sommes partielles en shared memory
  if (blockSize >= 512) {
	if (tib < 256) {
	  scumuls[ CFI(tib)].c1 += scumuls[ CFI(tib + 256) ].c1;
	  scumuls[ CFI(tib)].cx += scumuls[ CFI(tib + 256) ].cx;
	  scumuls[ CFI(tib)].cx2 += scumuls[CFI(tib + 256) ].cx2;
	}
	__syncthreads();
  }
 
  if (blockSize >= 256) {
	if (tib < 128) {
	  scumuls[ CFI(tib)].c1 += scumuls[ CFI(tib + 128) ].c1;
	  scumuls[ CFI(tib)].cx += scumuls[ CFI(tib + 128) ].cx;
	  scumuls[ CFI(tib)].cx2 += scumuls[CFI(tib + 128) ].cx2; 
	}
	__syncthreads();
  }
  if (blockSize >= 128) {
	if (tib <  64) {
	  scumuls[ CFI(tib)].c1 += scumuls[ CFI(tib +  64) ].c1;
	  scumuls[ CFI(tib)].cx += scumuls[ CFI(tib +  64) ].cx;
	  scumuls[CFI(tib)].cx2 += scumuls[CFI(tib +  64) ].cx2;	  
	}
	__syncthreads();
  }
  
  //32 threads <==> 1 warp
 volatile tcontribs * scontribs = scumuls ;
 if (tib < 32){
	 {
		scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib + 32) ].c1;
		scontribs[ CFI(tib)].cx += scontribs[ CFI(tib + 32) ].cx;
		scontribs[CFI(tib)].cx2 += scontribs[CFI(tib + 32) ].cx2;
	  }
 if (tib < 16)
	  {
		scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib + 16) ].c1;
		scontribs[ CFI(tib)].cx += scontribs[ CFI(tib + 16) ].cx;
		scontribs[CFI(tib)].cx2 += scontribs[CFI(tib + 16) ].cx2;
	  }
 if (tib < 8)
	  {
		scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  8) ].c1;
		scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  8) ].cx;
		scontribs[CFI(tib)].cx2 += scontribs[CFI(tib +  8) ].cx2;
	  }
 if (tib < 4)
   	{
	  scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  4) ].c1;
	  scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  4) ].cx;
	  scontribs[ CFI(tib)].cx2 +=scontribs[ CFI(tib +  4) ].cx2;
	}
 if (tib < 2)
   	{  
	  scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  2) ].c1;
	  scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  2) ].cx;
	  scontribs[CFI(tib)].cx2 += scontribs[CFI(tib +  2) ].cx2;
	 } 
 if (tib == 0)
   	{
	  scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  1) ].c1;
	  scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  1) ].cx;
	  scontribs[CFI(tib)].cx2 += scontribs[CFI(tib +  1) ].cx2;
  
  // resultat sommes partielles en gmem
  //if (tib == 0) {
	gsombloc[ blockIdx.x ] = scontribs[0].c1 ; 
	gsombloc[ blockIdx.x + gridDim.x ] = scontribs[0].cx ;
	gsombloc[ blockIdx.x + 2*gridDim.x ] = scontribs[0].cx2 ;
  	
 //code seg
	if (dy > 0 ) d_codes_x16[id_base_seg] = -1 ;
	if (dy < 0 ) d_codes_x16[id_base_seg] =  1 ;
	if (dy == 0) d_codes_x16[id_base_seg]=  0 ;
  }
 }
}

/*
 calcul des freeman et du centre de chaque segment de test
 a executer sur 'n_interval' blocs de 16 threads
 soit un thread par segment 
*/
__global__ void calcul_freemans_centre(uint2 * liste_pix, int * d_table_freeman, uint4 * d_freemans_x16){
	
  int id_segment = threadIdx.x ;    
  int id_freeman = ( blockDim.x*blockIdx.x + id_segment ) ;
  int id_base_pix = 5*id_freeman ;
  
  //calculs freemans, centre et code segment
  //1 uint4 par segment  
  int Dio, Djo, Dii, Dji;
      	
	//freeman out
	Dio = 1 + liste_pix[id_base_pix +1].x - liste_pix[id_base_pix].x ; 
	Djo = 1 + liste_pix[id_base_pix +1].y - liste_pix[id_base_pix].y ;
	d_freemans_x16[id_freeman].z = d_table_freeman[3*Dio + Djo] ;
	
	
	//freeman_in
	Dii = 1 + liste_pix[id_base_pix +4].x - liste_pix[id_base_pix +3].x ; 
	Dji = 1 + liste_pix[id_base_pix +4].y - liste_pix[id_base_pix +3].y ;
	d_freemans_x16[id_freeman].w = d_table_freeman[3*Dii + Dji] ;
	
	//centre 
	d_freemans_x16[id_freeman].x = liste_pix[id_base_pix +2].x ;
	d_freemans_x16[id_freeman].y = liste_pix[id_base_pix +2].y ;
	   
}


/*
  calcul des contribs 1, x et x2 des 16 segments
  autour du noeud
  a partir des sommes partielles des blocs
  1 bloc / 1 thread par segment 
*/

__global__ void somsom_full(uint64 * somblocs, int nb_nodes, unsigned int nb_bl_seg, uint64 * somsom){

  uint64 sdata[3];
  unsigned int seg = blockIdx.x ;
  unsigned int nb_seg = gridDim.x ;
  
  //un thread par segment
  {
	sdata[0] = 0;
	sdata[1] = 0;
	sdata[2] = 0;
  }

  for (int b=0; b < nb_bl_seg ; b++){ //voir atomicadd 64bits
	sdata[0] += somblocs[seg*nb_bl_seg + b];
	sdata[1] += somblocs[(seg + nb_seg)*nb_bl_seg + b];
	sdata[2] += somblocs[(seg + 2*nb_seg)*nb_bl_seg + b];
  }
  
  //totaux en gmem
  {
	somsom[3*seg] = sdata[0];
	somsom[3*seg + 1] = sdata[1];
	somsom[3*seg + 2] = sdata[2];
  }	  
}

/*
  version GPU de la fonction definie ds src/lib_math.c
*/
__device__ bool test_inf_gpu(double arg1, double arg2){
   if (arg2 > 0)
	 return (arg1 < (arg2*COEF_DECROI)) ;
  else
    return (arg1 < (arg2*INV_COEF_DECROI)) ;
}

/*
  version GPU de la fonction codage_niveau_gris_hyp_gaussienne
 */
__device__ double codage_gl_hyp_gauss(uint64 stat_sum_1, uint64 stat_sum_x, uint64 stat_sum_x2,
									  uint64 n_dim, uint64 SUM_X, uint64 SUM_X2){
  uint64 stat_sum_xe ;  /* somme des xn region exterieure */
  uint32 ne ;             /* nombre de pixel region exterieure */
  double sigi2, sige2; /* variance region interieure et exterieure */
  
  /* variance des valeurs des niveaux de gris a l'interieur du snake */
  sigi2 = 
    (double)stat_sum_x2/(double)stat_sum_1 - 
    ((double)stat_sum_x/stat_sum_1)*((double)stat_sum_x/stat_sum_1) ;

  /* variance des valeurs des niveaux de gris a l'exterieur du snake */
  ne = n_dim-stat_sum_1 ;
  stat_sum_xe = SUM_X - stat_sum_x ;
  sige2 =
    ((double)SUM_X2-stat_sum_x2)/ne - 
    ((double)stat_sum_xe/ne)*((double)stat_sum_xe/ne) ;
  
  if ((sigi2 > 0)&&(sige2 > 0))
  return  0.5*((double)stat_sum_1*log(sigi2) + (double)ne*log(sige2)) ;
  return -1 ;
}

/*
  soustrait, pour chaque intervalle [N1--Nx--N2] 
  les contributions des segments N1--Nx et Nx--N2
  
  A executer par nnodes blocs de 1 thread par intervalle
 */

__global__ void soustrait_aux_stats_2N_segments_noeud(snake_node_gpu * d_snake, int64 * d_stats_snake, int64 * d_stats_ref, 
													  t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2,
													  int * TABLE_CODAGE, uint32 l
													  )
{
  int nnodes = gridDim.x ;
  int id_nx, id_nprec, id_nprecprec, id_nsuiv ;
  int code_noeud, pos; 
  __shared__ int64 s_stats_snake[3] ;

  id_nx = blockIdx.x ;
  if (id_nx == 0) id_nprec = nnodes - 1;
  else id_nprec = id_nx - 1;
  if (id_nprec == 0) id_nprecprec = nnodes -1 ;
  else id_nprecprec = id_nprec - 1 ;
  if (id_nx == nnodes-1) id_nsuiv = 0;
  else id_nsuiv = id_nx + 1 ;
  
  //init avec les valeurs du contour actuel
  s_stats_snake[0] = d_stats_snake[0] ;
  s_stats_snake[1] = d_stats_snake[1] ;
  s_stats_snake[2] = d_stats_snake[2] ;

  /* segment Na -- Nx */
  if (d_snake[id_nprec].code_segment > 0)
    {
     s_stats_snake[0] -= d_snake[id_nprec].sum_1 ;
     s_stats_snake[1] -= d_snake[id_nprec].sum_x ;
     s_stats_snake[2] -= d_snake[id_nprec].sum_x2 ;
    }
  else if (d_snake[id_nprec].code_segment < 0)
    {
	 s_stats_snake[0] += d_snake[id_nprec].sum_1 ;
     s_stats_snake[1] += d_snake[id_nprec].sum_x ;
     s_stats_snake[2] += d_snake[id_nprec].sum_x2 ;
    }
  // else (code_segment_NaNx == 0), on ne fait rien
  
  /* gestion des pixels connectant les segments */
  /* pixel de depart du segment Na --> Nx */
  pos = d_snake[id_nprecprec].freeman_in*8 + d_snake[id_nprec].freeman_out ;
  code_noeud = TABLE_CODAGE[pos] ;
  pos = d_snake[id_nprec].posi*l + d_snake[id_nprec].posj ;
  if (code_noeud > 0)
    {
     s_stats_snake[0] -= 1 + d_snake[id_nprec].posj ;
     s_stats_snake[1] -= cumul_x[pos] ;
	 s_stats_snake[2] -= cumul_x2[pos];
    }
  else if (code_noeud < 0)
    {
	 s_stats_snake[0] += 1 + d_snake[id_nprec].posj ;
     s_stats_snake[1] += cumul_x[pos] ;
	 s_stats_snake[2] += cumul_x2[pos];
    }
  // else (code_noeud == 0), on ne fait rien
  
  /* -------------------------- */
  /* -------------------------- */
  /* segment Nx -- Nb */
  if ( d_snake[id_nx].code_segment > 0 )
    {
      // on soustrait la contribution 
	 s_stats_snake[0] -= d_snake[id_nx].sum_1 ;
     s_stats_snake[1] -= d_snake[id_nx].sum_x ;
     s_stats_snake[2] -= d_snake[id_nx].sum_x2 ;
    }
  else if ( d_snake[id_nx].code_segment < 0)
    {
	 s_stats_snake[0] += d_snake[id_nx].sum_1 ;
     s_stats_snake[1] += d_snake[id_nx].sum_x ;
     s_stats_snake[2] += d_snake[id_nx].sum_x2 ;
    }
  // else (code_segment_NxNb == 0), on ne fait rien
  
  /* gestion des pixels connectant les segments */
  /* pixel de depart du segment Nx --> Nb */
  pos = d_snake[id_nprec].freeman_in*8 + d_snake[id_nx].freeman_out ;
  code_noeud = TABLE_CODAGE[pos] ;
  pos = d_snake[id_nx].posi*l + d_snake[id_nx].posj ;
  if (code_noeud > 0)
    {
     s_stats_snake[0] -= 1 + d_snake[id_nx].posj  ;
     s_stats_snake[1] -= cumul_x[pos] ;
	 s_stats_snake[2] -= cumul_x2[pos];
    }
  else if (code_noeud < 0)
    {
	 s_stats_snake[0] += 1 + d_snake[id_nx].posj ;
     s_stats_snake[1] += cumul_x[pos] ;
	 s_stats_snake[2] += cumul_x2[pos];
    }
  // else (code_noeud == 0), on ne fait rien
  
  /* pixel d'arrivee du segment Nx --> Nb */
  pos = d_snake[id_nx].freeman_in*8 + d_snake[id_nsuiv].freeman_out ;
  code_noeud = TABLE_CODAGE[pos] ;
  pos = d_snake[id_nsuiv].posi*l +  d_snake[id_nsuiv].posj ;
  if (code_noeud > 0)
    {
	 s_stats_snake[0] -= 1 + d_snake[id_nsuiv].posj ;
     s_stats_snake[1] -= cumul_x[pos] ;
	 s_stats_snake[2] -= cumul_x2[pos];
    }
  else if (code_noeud < 0)
    {
	 s_stats_snake[0] += 1 + d_snake[id_nsuiv].posj ;
     s_stats_snake[1] += cumul_x[pos] ;
	 s_stats_snake[2] += cumul_x2[pos];
	 }
  // else (code_noeud == 0), on ne fait rien
  __syncthreads();
  
  d_stats_ref[3*id_nx] = s_stats_snake[0] ;
  d_stats_ref[3*id_nx + 1] = s_stats_snake[1] ;
  d_stats_ref[3*id_nx + 2] = s_stats_snake[2] ;
 
}


/*
  calcul des stats associees a chaque position de test

  EXEC : sur n_interval blocs de 8 threads
*/
__global__ void calcul_stats_full(snake_node_gpu * d_snake, int nnodes, bool pairs, int64 * d_stats_snake,
								  int64 * d_stats_ref, int64 * d_stats, uint64 * d_contribs,
								  uint4 * d_liste_points, int * code_segment, uint4 * d_freemans,
								  int * d_table_codes, t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2,
								  uint32 h, uint32 l, double * vrais, double * vrais_min, bool * move){

  int interval = blockIdx.x ;
  int seg = threadIdx.x ;
  int thread = seg;
  int id_nx, id_nprec, id_nprecprec, id_nsuiv, id_seg ;
  int code_noeud; 
  __shared__ int64 s_stats_snake[3*8] ;

  id_nx = 2*interval + !pairs ;
  if (id_nx == 0) id_nprec = nnodes - 1 ;
  else id_nprec = id_nx - 1 ;
  if (id_nprec == 0) id_nprecprec = nnodes -1 ;
  else id_nprecprec = id_nprec - 1 ;
  if (id_nx == nnodes-1) id_nsuiv = 0 ;
  else id_nsuiv = id_nx + 1 ;
  
  //chargement en smem , prevoir CFI car on a 24x64bits par blocs => conflits 
  s_stats_snake[3*thread + 0] = d_stats_ref[3*id_nx] ;
  s_stats_snake[3*thread + 1] = d_stats_ref[3*id_nx + 1] ;
  s_stats_snake[3*thread + 2] = d_stats_ref[3*id_nx + 2] ;
 
  
  //stats segments N1-Nx
  id_seg = 16*interval + seg ;
  if ( code_segment[id_seg] > 0 ){
	s_stats_snake[3*thread +0]  += d_contribs[3*id_seg] ;
	s_stats_snake[3*thread +1]  += d_contribs[3*id_seg + 1 ] ;
	s_stats_snake[3*thread +2]  += d_contribs[3*id_seg + 2 ] ;
  } else if ( code_segment[16*interval + seg] < 0 ) {
	s_stats_snake[3*thread +0]  -= d_contribs[3*id_seg] ;
	s_stats_snake[3*thread +1]  -= d_contribs[3*id_seg + 1] ;
	s_stats_snake[3*thread +2]  -= d_contribs[3*id_seg + 2] ;
  }
  
  //stats noeud N1(i1,j1)
  int fo_N1 = d_freemans[id_seg].z ;
  int fi_Nprecprec = d_snake[id_nprecprec].freeman_in ;
  int pos = d_snake[id_nprec].posi*l + d_snake[id_nprec].posj ;
   
  code_noeud = d_table_codes[fi_Nprecprec*8 + fo_N1];
  if (code_noeud > 0){
	s_stats_snake[3*thread +0] += 1 + d_snake[id_nprec].posj ;
	s_stats_snake[3*thread +1] += cumul_x[ pos ] ;
	s_stats_snake[3*thread +2] += cumul_x2[pos ] ;
  } else if (code_noeud < 0){
	s_stats_snake[3*thread +0] -= 1 + d_snake[id_nprec].posj ;
	s_stats_snake[3*thread +1] -= cumul_x[ pos ] ;
	s_stats_snake[3*thread +2] -= cumul_x2[pos ] ;
  }
  
  //stats noeud Nx
  int fo_Nx = d_freemans[id_seg + 8].z ;
  int fi_Nx = d_freemans[id_seg].w ;
  int Nxi = d_liste_points[8*id_nx + seg].x ;
  int Nxj = d_liste_points[8*id_nx + seg].y ;
  pos = Nxi*l + Nxj ;
  
  code_noeud = d_table_codes[fi_Nx*8 + fo_Nx];
  if (code_noeud > 0){
	s_stats_snake[3*thread +0] += 1 + Nxj ;
	s_stats_snake[3*thread +1] += cumul_x[ pos ] ;
	s_stats_snake[3*thread +2] += cumul_x2[pos ] ;
  }
  if (code_noeud < 0){
	s_stats_snake[3*thread +0] -= 1 + Nxj ;
	s_stats_snake[3*thread +1] -= cumul_x[ pos ] ;
	s_stats_snake[3*thread +2] -= cumul_x2[pos ] ;
  }
  
  //stats segments Nx-N2
  seg += 8;
  id_seg = 16*interval + seg ;
  if ( code_segment[id_seg] > 0 ){
	s_stats_snake[3*thread +0]  += d_contribs[3*id_seg] ;
	s_stats_snake[3*thread +1]  += d_contribs[3*id_seg + 1] ;
	s_stats_snake[3*thread +2]  += d_contribs[3*id_seg + 2] ;
  }
  if ( code_segment[id_seg] < 0 ) {
	s_stats_snake[3*thread +0]  -= d_contribs[3*id_seg] ;
	s_stats_snake[3*thread +1]  -= d_contribs[3*id_seg + 1] ;
	s_stats_snake[3*thread +2]  -= d_contribs[3*id_seg + 2] ;
  }
  
  //stats noeud N2(i2,j2)
  int fi_N2 = d_freemans[id_seg].w ;
  int fo_N2 = d_snake[id_nsuiv].freeman_out ;
  pos = d_snake[id_nsuiv].posi*l + d_snake[id_nsuiv].posj ;
  
  code_noeud = d_table_codes[fi_N2*8 + fo_N2];
  if (code_noeud > 0){
	s_stats_snake[3*thread +0] += 1 + d_snake[id_nsuiv].posj ;
	s_stats_snake[3*thread +1] += cumul_x[ pos ] ;
	s_stats_snake[3*thread +2] += cumul_x2[pos ] ;
  }
  if (code_noeud < 0){
	s_stats_snake[3*thread +0] -= 1 + d_snake[id_nsuiv].posj ;
	s_stats_snake[3*thread +1] -= cumul_x[ pos ] ;
	s_stats_snake[3*thread +2] -= cumul_x2[pos ] ;
  }
  
	  //TODO
	  //voir si on peut s'en passer
  d_stats[3*(8*interval + thread)] = s_stats_snake[3*thread +0]; 
  d_stats[3*(8*interval + thread) + 1] = s_stats_snake[3*thread +1];
  d_stats[3*(8*interval + thread) + 2] = s_stats_snake[3*thread +2];
  
  //codage hyp gaussienne  
  uint64 stat_sum_xe[8] ;    //somme des xn region exterieure 
  uint32 ne[8] ;             // nombre de pixels region exterieure 
  double sigi2[8], sige2[8]; //  carres des variances, regions interieure et exterieure 

  /* variance des valeurs des niveaux de gris a l'interieur du snake */
  sigi2[thread] =
	((double)s_stats_snake[3*thread +2]/(double)s_stats_snake[3*thread +0]) -
	((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]) ;
  
  /* variance des valeurs des niveaux de gris a l'exterieur du snake */
  ne[thread] = h*l-s_stats_snake[3*thread +0] ;
  stat_sum_xe[thread] = d_stats_snake[4] - s_stats_snake[3*thread +1] ;
  sige2[thread] =
	(double)(d_stats_snake[5]-s_stats_snake[3*thread +2])/(double)ne[thread] -
	((double)stat_sum_xe[thread]/ne[thread])*((double)stat_sum_xe[thread]/ne[thread]) ;

  if (sige2[thread]>0 && sigi2[thread]>0)
	vrais[8*interval + thread] = 0.5*((double)s_stats_snake[3*thread]*log(sigi2[thread]) + (double)ne[thread]*log(sige2[thread])) ;
  else
	vrais[8*interval + thread] = -1.0;

  if ( thread == 0 ){
	 //init move
	move[id_nx] = false ;
	int pos_optim = -1;
	double vrais_tmp = *vrais_min;
	for (int v=0; v < 8; v++){
	  if ( (vrais[8*interval + v] > 0) && (vrais[8*interval + v] < vrais_tmp*COEF_DECROI) ) {
		vrais_tmp = vrais[8*interval + v];
		pos_optim = v;
	  }
	}
	if (pos_optim >-1){
	 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) )
		{
		  /*maj data snake*/
		  move[id_nx] = true ;
		  //new position 
		  d_snake[id_nx].posi = d_liste_points[8*id_nx + pos_optim].x ;
		  d_snake[id_nx].posj = d_liste_points[8*id_nx + pos_optim].y ;
		  //nb pixels segment precedent
		  d_snake[id_nprec].nb_pixels = d_liste_points[8*id_nx + pos_optim].z ;
		  //nb pixels segment suivant
		  d_snake[id_nx].nb_pixels = d_liste_points[8*id_nx + pos_optim].w ;
		  //contribs segment precedent
		  d_snake[id_nprec].sum_1 = d_contribs[3*(16*interval + pos_optim)] ;
		  d_snake[id_nprec].sum_x = d_contribs[3*(16*interval + pos_optim) + 1] ;
		  d_snake[id_nprec].sum_x2 = d_contribs[3*(16*interval + pos_optim) + 2] ;
		  //contribs segment suivant
		  d_snake[id_nx].sum_1 = d_contribs[3*(16*interval + pos_optim + 8)] ;
		  d_snake[id_nx].sum_x = d_contribs[3*(16*interval + pos_optim + 8) + 1] ;
		  d_snake[id_nx].sum_x2 = d_contribs[3*(16*interval + pos_optim + 8) + 2] ;
		  //freemans segment precedent
		  d_snake[id_nprec].freeman_out = d_freemans[16*interval + pos_optim].z ;
		  d_snake[id_nprec].freeman_in = d_freemans[16*interval + pos_optim].w ;
		  //freemans segment suivant
		  d_snake[id_nx].freeman_out = d_freemans[16*interval + pos_optim + 8].z ;
		  d_snake[id_nx].freeman_in = d_freemans[16*interval + pos_optim + 8].w ;
		  //codes segment precedent
		  d_snake[id_nprec].code_segment = code_segment[16*interval + pos_optim] ;
		  //code segment suivant
		  d_snake[id_nx].code_segment = code_segment[16*interval + pos_optim + 8] ;
		  //centre segment precedent
		  d_snake[id_nprec].centre_i = d_freemans[16*interval + pos_optim ].x ;
		  d_snake[id_nprec].centre_j = d_freemans[16*interval + pos_optim ].y ;
		  //centre segment suivant 
		  d_snake[id_nx].centre_i = d_freemans[16*interval + pos_optim + 8].x ;
		  d_snake[id_nx].centre_j = d_freemans[16*interval + pos_optim + 8].y ;
		  
		}			   
	}
  }

  
}

__global__ void recalcul_stats_snake(snake_node_gpu * d_snake, int  nnodes, int64 * d_stats_snake, double * vrais_min,
									 t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2, int * TABLE_CODAGE, uint32 l
									 )
{
  
  int id_nx, id_nprec, id_nprecprec ;
  int code_noeud, code_segment, pos; 
  int64 s_stats_snake[3] ;
 
  //init stats en shared mem
  s_stats_snake[0] = 0 ;
  s_stats_snake[1] = 0 ;
  s_stats_snake[2] = 0 ;

    
  for (id_nx = 0; id_nx < nnodes; id_nx++)
    {
	  if (id_nx == 0) id_nprec = nnodes - 1;
	  else id_nprec = id_nx - 1;
	  if (id_nprec == 0) id_nprecprec = nnodes -1 ;
	  else id_nprecprec = id_nprec - 1 ;
      /* gestion des segments partant du noeud */
      /* vers le noeud suivant dans l'ordre trigo */
      code_segment = d_snake[id_nprec].code_segment ;
      if (code_segment > 0)
		{
		  /* on somme les contributions */
		  s_stats_snake[0] += d_snake[id_nprec].sum_1 ;
		  s_stats_snake[1] += d_snake[id_nprec].sum_x ;
		  s_stats_snake[2] += d_snake[id_nprec].sum_x2 ;
		}
      else if (code_segment < 0)
		{
		  /* on soustrait les contributions */
		  s_stats_snake[0] -= d_snake[id_nprec].sum_1 ;
		  s_stats_snake[1] -= d_snake[id_nprec].sum_x ;
		  s_stats_snake[2] -= d_snake[id_nprec].sum_x2 ;
		}
      // else (code_segment == 0), on ne fait rien
      /* gestion des pixels connectant les segments */
	  pos = d_snake[id_nprecprec].freeman_in*8 + d_snake[id_nprec].freeman_out ;
      code_noeud = TABLE_CODAGE[pos] ;
	  pos = d_snake[id_nprec].posi*l + d_snake[id_nprec].posj ;
   
      if (code_noeud > 0)
		{
		  /* on somme les contributions */
		  s_stats_snake[0] += 1 + d_snake[id_nprec].posj ;
		  s_stats_snake[1] += cumul_x[pos] ;
		  s_stats_snake[2] += cumul_x2[pos] ;
		}
      else if (code_noeud < 0)
		{
		  /* on soustrait les contributions */
		  s_stats_snake[0] -= 1 + d_snake[id_nprec].posj ;
		  s_stats_snake[1] -= cumul_x[pos] ;
		  s_stats_snake[2] -= cumul_x2[pos] ;
		}
      // else (code_pixel == 0), on ne fait rien
    }
  d_stats_snake[0] = s_stats_snake[0] ;
  d_stats_snake[1] = s_stats_snake[1] ;
  d_stats_snake[2] = s_stats_snake[2] ;
  
  *vrais_min = codage_gl_hyp_gauss(s_stats_snake[0], s_stats_snake[1], s_stats_snake[2],
								   d_stats_snake[3], d_stats_snake[4], d_stats_snake[5]);
}


__global__ void ajoute_noeuds(snake_node_gpu * snake, snake_node_gpu * snake_tmp, int nnodes, int seuil, int * new_nb_nodes){

  volatile snake_node_gpu * st = snake_tmp ;

 int id_cpy = 0;  
 for (int id_nx=0; id_nx < nnodes; id_nx++){
	//position du noeud existant
	st[id_cpy].posi =  snake[id_nx].posi ;
    st[id_cpy].posj =  snake[id_nx].posj ;
        
	id_cpy++ ;
	
	if ( snake[id_nx].nb_pixels > seuil)
	{
		//position du nouveau noeud
		st[id_cpy].posi = snake[id_nx].centre_i ;
		st[id_cpy].posj = snake[id_nx].centre_j ;
		id_cpy++ ;
	}	
 }
 for( int node=0; node<id_cpy; node++ ){
   snake[node].posi = st[node].posi ;
   snake[node].posj = st[node].posj ;
 }
*new_nb_nodes = id_cpy-nnodes ;
}
  




__global__ void recalcul_contribs_segments_snake(snake_node_gpu * d_snake, int nb_nodes,
												 t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2, 
												 int l, uint2 * liste_pix, uint64 * gsombloc )
{
  // indices des elements
  int blockSize = blockDim.x ;
  int tib = threadIdx.x ;
  int nblocs_seg =  gridDim.x / nb_nodes ;
  int idx = blockDim.x*blockIdx.x + threadIdx.x ;
  int segment = blockIdx.x / nblocs_seg ;
  int tis = idx - segment*nblocs_seg*blockDim.x ;

  //tab pour contribs pixels 
  extern __shared__ tcontribs scumuls[];

  //indices des noeuds
  uint x1, y1, x2, y2 ;
  int n1, n2 ;
  uint2 p ;
  int xsuiv, xprec ;
  
  n1 = segment ;
  n2 = segment +1 ;
  //gestion du bouclage du snake
  if (n2 >= nb_nodes) n2 = 0 ;
  
  //affectation des differentes positions aux différents segments 'blocs de threads'
  	x1 = d_snake[n1].posj ;
	y1 = d_snake[n1].posi ;
	x2 = d_snake[n2].posj ;
	y2 = d_snake[n2].posi ;
  
  //params des deplacements
  int dx=x2-x1;
  int dy=y2-y1;
  uint abs_dx = ABS(dx);
  uint abs_dy = ABS(dy);
  uint nb_pix = abs_dy>abs_dx?(abs_dy+1):(abs_dx+1); // alternative -> lecture ds liste_points[]
  int incx=0, incy=0;
  
  
  //calcul liste des pixels du segment (x1,y1)-(x2,y2)  
  if (dy > 0) incy=1; else incy=-1 ;
  if (dx > 0) incx=1; else incx=-1 ;
  
  if (tis < nb_pix){
	if (abs_dy > abs_dx){
	  //1 thread par ligne
	  double k = (double)dx/dy ;
	  p.x = y1 + incy*tis ;
	  p.y = x1 + floor((double)incy*k*tis+0.5) ;

	} else {
	  //1 thread par colonne	  
	  double k=(double)dy/dx ;
	  p.x = y1 + floor((double)(incx*k*tis)+0.5) ;
	  p.y =  x1 + incx*tis ;
	  if ( tis > 0 ){ 
		xsuiv = y1 + floor((double)(incx*k*(tis+1))+0.5) ;
		xprec = y1 + floor((double)(incx*k*(tis-1))+0.5) ;
	  }

	}
	if (tis == 0)  liste_pix[5*segment] = make_uint2(p.x, p.y) ;
	if (tis == 1)  liste_pix[5*segment +1] = make_uint2(p.x,p.y) ;
    if (tis == nb_pix/2) liste_pix[5*segment +2] = make_uint2(p.x,p.y) ;
    if (tis == nb_pix-2) liste_pix[5*segment +3] = make_uint2(p.x,p.y) ;
    if (tis == nb_pix-1) liste_pix[5*segment +4] = make_uint2(p.x,p.y) ;
 
  }
  __syncthreads();
    
  //calcul contribs individuelles des pixels
  
  if ( (tis >0) && (tis < nb_pix-1)
	   && ( ( (abs_dy <= abs_dx) && ( ( xprec > p.x) || ( xsuiv > p.x)) )
			|| (abs_dy > abs_dx) ) )
	{
	int pos = p.x * l + p.y ;
	scumuls[ CFI(tib)].c1 = 1+p.y;        
	scumuls[ CFI(tib)].cx = cumul_x[ pos ] ;
	scumuls[CFI(tib)].cx2 = cumul_x2[ pos ];
  } else {
	scumuls[ CFI(tib)].c1 = 0;
	scumuls[ CFI(tib)].cx = 0;
	scumuls[CFI(tib)].cx2 = 0;
  }
  
   __syncthreads();
  //somme des contribs individuelles 
  // unroll des  sommes partielles en smem
  
  if (blockSize >= 512) {
	if (tib < 256) {
	  scumuls[ CFI(tib)].c1 += scumuls[ CFI(tib + 256) ].c1;
	  scumuls[ CFI(tib)].cx += scumuls[ CFI(tib + 256) ].cx;
	  scumuls[CFI(tib)].cx2 += scumuls[CFI(tib + 256) ].cx2;
	}
	__syncthreads();
  }
 
  if (blockSize >= 256) {
	if (tib < 128) {
	  scumuls[ CFI(tib)].c1 += scumuls[ CFI(tib + 128) ].c1;
	  scumuls[ CFI(tib)].cx += scumuls[ CFI(tib + 128) ].cx;
	  scumuls[CFI(tib)].cx2 += scumuls[CFI(tib + 128) ].cx2; 
	}
	__syncthreads();
  }
  if (blockSize >= 128) {
	if (tib <  64) {
	  scumuls[ CFI(tib)].c1 += scumuls[ CFI(tib +  64) ].c1;
	  scumuls[ CFI(tib)].cx += scumuls[ CFI(tib +  64) ].cx;
	  scumuls[CFI(tib)].cx2 += scumuls[ CFI(tib +  64) ].cx2;	  
	}
	__syncthreads();
  }
  
  //32 threads <==> 1 warp
  volatile tcontribs * scontribs = scumuls ;
  if (tib < 32)
	{
	  {
		scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib + 32) ].c1;
		scontribs[ CFI(tib)].cx += scontribs[ CFI(tib + 32) ].cx;
		scontribs[ CFI(tib)].cx2 +=scontribs[ CFI(tib + 32) ].cx2;
	  }
	  if (tib<16)
	  {
		scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib + 16) ].c1;
		scontribs[ CFI(tib)].cx += scontribs[ CFI(tib + 16) ].cx;
		scontribs[ CFI(tib)].cx2 +=scontribs[ CFI(tib + 16) ].cx2;
	  }
	  if (tib<8)
	  {
		scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  8) ].c1;
		scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  8) ].cx;
		scontribs[ CFI(tib)].cx2 +=scontribs[ CFI(tib +  8) ].cx2;}
	if (tib<4){	
	  scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  4) ].c1;
	  scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  4) ].cx;
	  scontribs[ CFI(tib)].cx2 +=scontribs[ CFI(tib +  4) ].cx2;
	  }
	  if (tib<2){
	  scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  2) ].c1;
	  scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  2) ].cx;
	  scontribs[ CFI(tib)].cx2 +=scontribs[ CFI(tib +  2) ].cx2;
	  }
	  if (tib==0){
	  scontribs[ CFI(tib)].c1 += scontribs[ CFI(tib +  1) ].c1;
	  scontribs[ CFI(tib)].cx += scontribs[ CFI(tib +  1) ].cx;
	  scontribs[ CFI(tib)].cx2 +=scontribs[ CFI(tib +  1) ].cx2;
  // resultat sommes partielles en gmem
  //if (tib == 0) {
	gsombloc[ blockIdx.x ] = scontribs[0].c1; 
	gsombloc[ blockIdx.x + gridDim.x ] = scontribs[0].cx;
	gsombloc[ blockIdx.x + 2*gridDim.x ] = scontribs[0].cx2;

  //calculs code segment
	//code seg
	if (dy > 0 ) d_snake[segment].code_segment = -1 ;
	if (dy < 0 ) d_snake[segment].code_segment = 1 ;
	if (dy == 0) d_snake[segment].code_segment = 0 ;
  }
	}
}

__global__ void recalcul_freemans_centre(snake_node_gpu * snake, uint2 * liste_pix, int * d_table_freeman){
	
  int id_segment = blockIdx.x ;    
  int id_base_pix = 5*id_segment ;    

  //calculs freemans, centre et code segment
  //1 uint4 par segment  
  int Dio, Djo, Dii, Dji;
      	
	//freeman out
	Dio = 1 + liste_pix[id_base_pix +1].x - liste_pix[id_base_pix].x ; 
	Djo = 1 + liste_pix[id_base_pix +1].y - liste_pix[id_base_pix].y ;
	snake[id_segment].freeman_out = d_table_freeman[3*Dio + Djo] ;
		
	//freeman_in
	Dii = 1 + liste_pix[id_base_pix +4].x - liste_pix[id_base_pix +3].x ; 
	Dji = 1 + liste_pix[id_base_pix +4].y - liste_pix[id_base_pix +3].y ;
	snake[id_segment].freeman_in = d_table_freeman[3*Dii + Dji] ;
	
	//centre 
	snake[id_segment].centre_i = liste_pix[id_base_pix +2].x ;
	snake[id_segment].centre_j = liste_pix[id_base_pix +2].y ;
	   
}


/*
  sommme des contribs par bloc -> contribs segment, pour le snake

  execution sur : 1bloc / 1 thread par segment
 */

__global__ void resomsom_snake(uint64 * somblocs, int nb_nodes, unsigned int nb_bl_seg, snake_node_gpu * d_snake){

  uint64 sdata[3];
  unsigned int seg = blockIdx.x ;
  
  //un thread par segment
  {
	sdata[0] = 0;
	sdata[1] = 0;
	sdata[2] = 0;
  }

  for (int b=0; b < nb_bl_seg ; b++){
	sdata[0] += somblocs[seg*nb_bl_seg + b];
	sdata[1] += somblocs[(seg + nb_nodes)*nb_bl_seg + b];
	sdata[2] += somblocs[(seg + 2*nb_nodes)*nb_bl_seg + b];
  }
  
  //totaux en gmem
  {
	d_snake[seg].sum_1 = sdata[0];
	d_snake[seg].sum_x = sdata[1];
	d_snake[seg].sum_x2 = sdata[2];
  }	  
}