#include "lib_kernels_maths.h"

__device__ void calcul_indices_prec_suiv(int nb_nodes, int nx, int& nprec, int& nsuiv){
  nprec = (nx > 0)? (nx - 1):(nb_nodes - 1);
  nsuiv = (nx < nb_nodes - 1)? (nx + 1):(0);
}

__device__ double codage_niveau_gris_hyp_gaussienne_gpu(uint32 stat_sum_1, uint64 stat_sum_x,
														uint64 stat_sum_x2, uint32 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_x)/(double)((uint64)stat_sum_1*(uint64)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)/(double)ne - 
    (double)(stat_sum_xe*stat_sum_xe)/(double)((uint64)ne*(uint64)ne) ;
  
  if ((sigi2 <= 0)|(sige2 <= 0)) return -1;
    
  return  0.5*((double)stat_sum_1*log(sigi2) + (double)ne*log(sige2)) ;
}

__device__ double codage_gl_hyp_gaussienne(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_x)/(double)((uint64)stat_sum_1*(uint64)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)/(double)ne - 
    (double)(stat_sum_xe*stat_sum_xe)/(double)((uint64)ne*(uint64)ne) ;
  
  if ((sigi2 > 0)|(sige2 > 0))
  return  0.5*((double)stat_sum_1*log(sigi2) + (double)ne*log(sige2)) ;
  return -1 ;
}

__device__ inline unsigned int nextPow2_gpu( unsigned int x ) {
  --x;
  x |= x >> 1;
  x |= x >> 2;
  x |= x >> 4;
  x |= x >> 8;
  x |= x >> 16;
  return ++x;
}


__device__ inline int sign_diff_strict(int val1, int val2)
{
  if (val1 > 0) 
    {
      if (val2 >= 0) return 0 ;
      else return 1 ;
    }
  else 
    if (val1 < 0)
      {
	if (val2 <= 0) return 0 ;
	else  return 1 ;
      }
    else
      return 0 ;/* val1 == 0 */
}

__device__ inline int sign_diff_ou_egal_zero(int val1, int val2)
{
  if (val1 > 0) 
    {
      if (val2 > 0) return 0 ;
      else return 1 ;
    }
  else 
    if (val1 < 0)
      {
		if (val2 < 0) return 0 ;
		else  return 1 ;
      }
    else
      return 1 ;/* val1 == 0 */
}


__device__ inline int sinus_triangle(int Ai, int Aj, int Bi, int Bj, int Ci, int Cj)
{
  return (((Bi-Ai)*(Cj-Aj)) - ((Ci-Ai)*(Bj-Aj))) ;
}

__device__ inline int test_croisement_large(uint32 Ai, uint32 Aj, uint32 Bi, uint32 Bj,
											uint32 Ui, uint32 Uj, uint32 Vi, uint32 Vj)
{
  if (sign_diff_strict(sinus_triangle(Ai, Aj, Bi, Bj, Ui, Uj),     // Bi-Ai, Uj-Aj, Ui-Ai, Bj-Aj
		       sinus_triangle(Ai, Aj, Bi, Bj, Vi, Vj)))    // Bi-Ai, Vj-Aj, Vi-Ai, Bj-Aj
    if (sign_diff_strict(sinus_triangle(Ui, Uj, Vi, Vj, Ai, Aj),   // Vi-Ui, Aj-Uj, Ai-Ui, Vj-Uj
			 sinus_triangle(Ui, Uj, Vi, Vj, Bi, Bj)))  // Vi-Ui, Bj-Uj, Bi-Ui, Vj-Uj
       return 1 ;
  return 0 ;
}

__device__ inline int test_croisement_strict(uint32 Ai, uint32 Aj, uint32 Bi, uint32 Bj,
											 uint32 Ui, uint32 Uj, uint32 Vi, uint32 Vj )
{
  if (sign_diff_ou_egal_zero(sinus_triangle(Ai, Aj, Bi, Bj, Ui, Uj),
							 sinus_triangle(Ai, Aj, Bi, Bj, Vi, Vj)))  
    if (sign_diff_ou_egal_zero(sinus_triangle(Ui, Uj, Vi, Vj, Ai, Aj),
							   sinus_triangle(Ui, Uj, Vi, Vj, Bi, Bj)))
      return 1 ;
  return 0 ;
}

__global__ void kernel_test_croisement_move_seg_strict(struct snake_node_gpu *d_snake, uint32 nx, uint32 Nxi, uint32 Nxj, int Nb_nodes, bool * croist)
{
  int idx = blockDim.x*blockIdx.x + threadIdx.x, idx_prec, idx_suiv ;
  int na, nb, naa, nbb  ;
  int Nai, Naj, Nbi, Nbj ;
  bool croistil = false ;
  
  calcul_indices_prec_suiv(Nb_nodes, nx, na, nb);
  calcul_indices_prec_suiv(Nb_nodes, na, naa, nbb);

  //coord. nodes extremites des segments a tester
  Nai = d_snake[na].posi ;
  Naj = d_snake[na].posj ;
  Nbi = d_snake[nb].posi ;
  Nbj = d_snake[nb].posj ;
  
  if (idx == nb) //premier segment apres
	{
	  calcul_indices_prec_suiv(Nb_nodes, idx, idx_prec, idx_suiv);
	  if ( test_croisement_large(Nxi, Nxj, Nbi, Nbj,
								  Nbi, Nbj, d_snake[idx_suiv].posi, d_snake[idx_suiv].posj))
		croistil = true ;
	  else if ( test_croisement_strict(Nai, Naj, Nxi, Nxj,
									   Nbi, Nbj, d_snake[idx_suiv].posi, d_snake[idx_suiv].posj) )
		croistil = true ;
	}
  else if ( idx == naa ) //premier segment avant
	{
	  if ( test_croisement_large(d_snake[naa].posi, d_snake[naa].posj, Nai, Naj,
								 Nai, Naj, Nxi, Nxj) )
		croistil = true ;
	  else if ( test_croisement_strict(d_snake[naa].posi, d_snake[naa].posj, Nai, Naj,
									   Nxi, Nxj, Nbi, Nbj))
		croistil = true ;
	}
  else if ( (idx != nx) && (idx != na) ) //les autres segments
	{
	  calcul_indices_prec_suiv(Nb_nodes, idx, idx_prec, idx_suiv);
	  if ( test_croisement_strict(Nai, Naj, Nxi, Nxj,
								  d_snake[idx].posi, d_snake[idx].posj,
								  d_snake[idx_suiv].posi, d_snake[idx_suiv].posj) )
		croistil = true ;
	  else if ( test_croisement_strict(Nxi, Nxj, Nbi, Nbj,
									   d_snake[idx].posi, d_snake[idx].posj,
									   d_snake[idx_suiv].posi, d_snake[idx_suiv].posj) )
		croistil = true ;
	}
	 
  croist[idx] = croistil ;
}

__device__ bool test_croisement(struct snake_node_gpu *d_snake, uint32 nx, uint32 Nxi, uint32 Nxj, int Nb_nodes)
{
  int idx_prec, idx_suiv ;
  int na, nb, naa, nbb  ;
  int Nai, Naj, Nbi, Nbj ;
  bool croistil = false ;

  //coord. nodes extremites des segments a tester
  calcul_indices_prec_suiv(Nb_nodes, nx, na, nb);
  calcul_indices_prec_suiv(Nb_nodes, na, naa, nbb);

  Nai = d_snake[na].posi ;
  Naj = d_snake[na].posj ;
  Nbi = d_snake[nb].posi ;
  Nbj = d_snake[nb].posj ;

  //parcours du snake
  for (int idx=0; idx < Nb_nodes; idx++){
	if (idx == nb) //premier segment apres
	  {
		calcul_indices_prec_suiv(Nb_nodes, idx, idx_prec, idx_suiv);
		if ( test_croisement_large(Nxi, Nxj, Nbi, Nbj,
								   Nbi, Nbj, d_snake[idx_suiv].posi, d_snake[idx_suiv].posj))
		  croistil = true ;
		else if ( test_croisement_strict(Nai, Naj, Nxi, Nxj,
									   Nbi, Nbj, d_snake[idx_suiv].posi, d_snake[idx_suiv].posj) )
		  croistil = true ;
	  }
	else if ( idx == naa ) //premier segment avant
	  {
		if ( test_croisement_large(d_snake[naa].posi, d_snake[naa].posj, Nai, Naj,
								   Nai, Naj, Nxi, Nxj) )
		  croistil = true ;
		else if ( test_croisement_strict(d_snake[naa].posi, d_snake[naa].posj, Nai, Naj,
										 Nxi, Nxj, Nbi, Nbj))
		  croistil = true ;
	  }
	else if ( (idx != nx) && (idx != na) ) //les autres segments
	  {
		calcul_indices_prec_suiv(Nb_nodes, idx, idx_prec, idx_suiv);
		if ( test_croisement_strict(Nai, Naj, Nxi, Nxj,
									d_snake[idx].posi, d_snake[idx].posj,
									d_snake[idx_suiv].posi, d_snake[idx_suiv].posj) )
		  croistil = true ;
		else if ( test_croisement_strict(Nxi, Nxj, Nbi, Nbj,
										 d_snake[idx].posi, d_snake[idx].posj,
										 d_snake[idx_suiv].posi, d_snake[idx_suiv].posj) )
		  croistil = true ;
	  }
  }	 
  return ( croistil );
}