src/lib_kernels_maths.cu

   1
   2 __device__ void calcul_indices_prec_suiv(int nb_nodes, int nx, int& nprec, int& nsuiv){
   3   nprec = (nx > 0)? (nx - 1):(nb_nodes - 1);
   4   nsuiv = (nx < nb_nodes - 1)? (nx + 1):(0);
   5 }
   6
   7 __device__ double codage_niveau_gris_hyp_gaussienne_gpu(uint32 stat_sum_1, uint64 stat_sum_x,
   8                                                                                                                 uint64 stat_sum_x2, uint32 n_dim,
   9                                                                                                                 uint64 SUM_X, uint64 SUM_X2)
  10 {
  11   uint64 stat_sum_xe ;  /* somme des xn region exterieure */
  12   uint32 ne ;             /* nombre de pixel region exterieure */
  13   double sigi2, sige2; /* variance region interieure et exterieure */
  14
  15   /* variance des valeurs des niveaux de gris a l'interieur du snake */
  16   sigi2 =
  17     (double)stat_sum_x2/(double)stat_sum_1 -
  18     (double)(stat_sum_x*stat_sum_x)/(double)((uint64)stat_sum_1*(uint64)stat_sum_1) ;
  19
  20   /* variance des valeurs des niveaux de gris a l'exterieur du snake */
  21   ne = n_dim-stat_sum_1 ;
  22   stat_sum_xe = SUM_X - stat_sum_x ;
  23   sige2 =
  24     (double)(SUM_X2-stat_sum_x2)/(double)ne -
  25     (double)(stat_sum_xe*stat_sum_xe)/(double)((uint64)ne*(uint64)ne) ;
  26
  27   if ((sigi2 <= 0)|(sige2 <= 0)) return -1;
  28
  29   return  0.5*((double)stat_sum_1*log(sigi2) + (double)ne*log(sige2)) ;
  30 }
  31
  32 __device__ double codage_gl_hyp_gaussienne(uint64 stat_sum_1, uint64 stat_sum_x, uint64 stat_sum_x2,
  33                                                                                    uint64 n_dim, uint64 SUM_X, uint64 SUM_X2){
  34   uint64 stat_sum_xe ;  /* somme des xn region exterieure */
  35   uint32 ne ;             /* nombre de pixel region exterieure */
  36   double sigi2, sige2; /* variance region interieure et exterieure */
  37
  38   /* variance des valeurs des niveaux de gris a l'interieur du snake */
  39   sigi2 =
  40     (double)stat_sum_x2/(double)stat_sum_1 -
  41     (double)(stat_sum_x*stat_sum_x)/(double)((uint64)stat_sum_1*(uint64)stat_sum_1) ;
  42
  43   /* variance des valeurs des niveaux de gris a l'exterieur du snake */
  44   ne = n_dim-stat_sum_1 ;
  45   stat_sum_xe = SUM_X - stat_sum_x ;
  46   sige2 =
  47     (double)(SUM_X2-stat_sum_x2)/(double)ne -
  48     (double)(stat_sum_xe*stat_sum_xe)/(double)((uint64)ne*(uint64)ne) ;
  49
  50   if ((sigi2 > 0)|(sige2 > 0))
  51   return  0.5*((double)stat_sum_1*log(sigi2) + (double)ne*log(sige2)) ;
  52   return -1 ;
  53 }
  54
  55 __device__ inline unsigned int nextPow2_gpu( unsigned int x ) {
  56   --x;
  57   x |= x >> 1;
  58   x |= x >> 2;
  59   x |= x >> 4;
  60   x |= x >> 8;
  61   x |= x >> 16;
  62   return ++x;
  63 }
  64
  65
  66 __device__ inline int sign_diff_strict(int val1, int val2)
  67 {
  68   if (val1 > 0)
  69     {
  70       if (val2 >= 0) return 0 ;
  71       else return 1 ;
  72     }
  73   else
  74     if (val1 < 0)
  75       {
  76         if (val2 <= 0) return 0 ;
  77         else  return 1 ;
  78       }
  79     else
  80       return 0 ;/* val1 == 0 */
  81 }
  82
  83 __device__ inline int sign_diff_ou_egal_zero(int val1, int val2)
  84 {
  85   if (val1 > 0)
  86     {
  87       if (val2 > 0) return 0 ;
  88       else return 1 ;
  89     }
  90   else
  91     if (val1 < 0)
  92       {
  93                 if (val2 < 0) return 0 ;
  94                 else  return 1 ;
  95       }
  96     else
  97       return 1 ;/* val1 == 0 */
  98 }
  99
 100
 101 __device__  int sinus_triangle(int Ai, int Aj, int Bi, int Bj, int Ci, int Cj)
 102 {
 103   return (((Bi-Ai)*(Cj-Aj)) - ((Ci-Ai)*(Bj-Aj))) ;
 104 }
 105
 106 __device__  int test_croisement_large(uint32 Ai, uint32 Aj, uint32 Bi, uint32 Bj,
 107                                                                                         uint32 Ui, uint32 Uj, uint32 Vi, uint32 Vj)
 108 {
 109   if (sign_diff_strict(sinus_triangle(Ai, Aj, Bi, Bj, Ui, Uj),     // Bi-Ai, Uj-Aj, Ui-Ai, Bj-Aj
 110                        sinus_triangle(Ai, Aj, Bi, Bj, Vi, Vj)))    // Bi-Ai, Vj-Aj, Vi-Ai, Bj-Aj
 111     if (sign_diff_strict(sinus_triangle(Ui, Uj, Vi, Vj, Ai, Aj),   // Vi-Ui, Aj-Uj, Ai-Ui, Vj-Uj
 112                          sinus_triangle(Ui, Uj, Vi, Vj, Bi, Bj)))  // Vi-Ui, Bj-Uj, Bi-Ui, Vj-Uj
 113        return 1 ;
 114   return 0 ;
 115 }
 116
 117 __device__  int test_croisement_strict(uint32 Ai, uint32 Aj, uint32 Bi, uint32 Bj,
 118                                                                                          uint32 Ui, uint32 Uj, uint32 Vi, uint32 Vj )
 119 {
 120   if (sign_diff_ou_egal_zero(sinus_triangle(Ai, Aj, Bi, Bj, Ui, Uj),
 121                                                          sinus_triangle(Ai, Aj, Bi, Bj, Vi, Vj)))
 122     if (sign_diff_ou_egal_zero(sinus_triangle(Ui, Uj, Vi, Vj, Ai, Aj),
 123                                                            sinus_triangle(Ui, Uj, Vi, Vj, Bi, Bj)))
 124       return 1 ;
 125   return 0 ;
 126 }
 127
 128 __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)
 129 {
 130   int idx = blockDim.x*blockIdx.x + threadIdx.x, idx_prec, idx_suiv ;
 131   int na, nb, naa, nbb  ;
 132   int Nai, Naj, Nbi, Nbj ;
 133   bool croistil = false ;
 134
 135   calcul_indices_prec_suiv(Nb_nodes, nx, na, nb);
 136   calcul_indices_prec_suiv(Nb_nodes, na, naa, nbb);
 137
 138   //coord. nodes extremites des segments a tester
 139   Nai = d_snake[na].posi ;
 140   Naj = d_snake[na].posj ;
 141   Nbi = d_snake[nb].posi ;
 142   Nbj = d_snake[nb].posj ;
 143
 144   if (idx == nb) //premier segment apres
 145         {
 146           calcul_indices_prec_suiv(Nb_nodes, idx, idx_prec, idx_suiv);
 147           if ( test_croisement_large(Nxi, Nxj, Nbi, Nbj,
 148                                                                   Nbi, Nbj, d_snake[idx_suiv].posi, d_snake[idx_suiv].posj))
 149                 croistil = true ;
 150           else if ( test_croisement_strict(Nai, Naj, Nxi, Nxj,
 151                                                                            Nbi, Nbj, d_snake[idx_suiv].posi, d_snake[idx_suiv].posj) )
 152                 croistil = true ;
 153         }
 154   else if ( idx == naa ) //premier segment avant
 155         {
 156           if ( test_croisement_large(d_snake[naa].posi, d_snake[naa].posj, Nai, Naj,
 157                                                                  Nai, Naj, Nxi, Nxj) )
 158                 croistil = true ;
 159           else if ( test_croisement_strict(d_snake[naa].posi, d_snake[naa].posj, Nai, Naj,
 160                                                                            Nxi, Nxj, Nbi, Nbj))
 161                 croistil = true ;
 162         }
 163   else if ( (idx != nx) && (idx != na) ) //les autres segments
 164         {
 165           calcul_indices_prec_suiv(Nb_nodes, idx, idx_prec, idx_suiv);
 166           if ( test_croisement_strict(Nai, Naj, Nxi, Nxj,
 167                                                                   d_snake[idx].posi, d_snake[idx].posj,
 168                                                                   d_snake[idx_suiv].posi, d_snake[idx_suiv].posj) )
 169                 croistil = true ;
 170           else if ( test_croisement_strict(Nxi, Nxj, Nbi, Nbj,
 171                                                                            d_snake[idx].posi, d_snake[idx].posj,
 172                                                                            d_snake[idx_suiv].posi, d_snake[idx_suiv].posj) )
 173                 croistil = true ;
 174         }
 175
 176   croist[idx] = croistil ;
 177 }
 178
 179 __device__ bool croisement(snake_node_gpu *d_snake, uint32 nx, uint32 Nxi, uint32 Nxj, int Nb_nodes)
 180 {
 181   int idx_prec, idx_suiv ;
 182   int na, nb, naa, nbb  ;
 183   int Nai, Naj, Nbi, Nbj ;
 184   bool croistil = false ;
 185
 186   //coord. nodes extremites des segments a tester
 187   calcul_indices_prec_suiv(Nb_nodes, nx, na, nb);
 188   calcul_indices_prec_suiv(Nb_nodes, na, naa, nbb);
 189
 190   Nai = d_snake[na].posi ;
 191   Naj = d_snake[na].posj ;
 192   Nbi = d_snake[nb].posi ;
 193   Nbj = d_snake[nb].posj ;
 194
 195   //parcours du snake
 196   for (int idx=0; idx < Nb_nodes; idx++){
 197         if (idx == nb) //premier segment apres
 198           {
 199                 calcul_indices_prec_suiv(Nb_nodes, idx, idx_prec, idx_suiv);
 200                 if ( test_croisement_large(Nxi, Nxj, Nbi, Nbj, Nbi, Nbj, d_snake[idx_suiv].posi, d_snake[idx_suiv].posj) )
 201                   {
 202                         croistil = true ;
 203                         break ;
 204                   }
 205                 if ( test_croisement_strict(Nai, Naj, Nxi, Nxj, Nbi, Nbj, d_snake[idx_suiv].posi, d_snake[idx_suiv].posj) )
 206                   {
 207                         croistil = true ;
 208                         break ;
 209                   }
 210           }
 211         if ( idx == naa ) //premier segment avant
 212           {
 213                 if ( test_croisement_large(d_snake[naa].posi, d_snake[naa].posj, Nai, Naj, Nai, Naj, Nxi, Nxj) )
 214                   {
 215                         croistil = true ;
 216                         break ;
 217                   }
 218                 if ( test_croisement_strict(d_snake[naa].posi, d_snake[naa].posj, Nai, Naj, Nxi, Nxj, Nbi, Nbj))
 219                   {
 220                         croistil = true ;
 221                         break ;
 222                   }
 223           }
 224         if ( (idx != nx) && (idx != na) && (idx != nb) && (idx != naa) ) //les autres segments
 225           {
 226                 calcul_indices_prec_suiv(Nb_nodes, idx, idx_prec, idx_suiv);
 227                 if ( test_croisement_strict(Nai, Naj, Nxi, Nxj,
 228                                                                         d_snake[idx].posi, d_snake[idx].posj,
 229                                                                         d_snake[idx_suiv].posi, d_snake[idx_suiv].posj) )
 230                   {
 231                         croistil = true ;
 232                         break ;
 233                   }
 234                 if ( test_croisement_strict(Nxi, Nxj, Nbi, Nbj,
 235                                                                         d_snake[idx].posi, d_snake[idx].posj,
 236                                                                         d_snake[idx_suiv].posi, d_snake[idx_suiv].posj) )
 237                   {
 238                         croistil = true ;
 239                         break ;
 240                   }
 241           }
 242   }
 243   return croistil;
 244 }