src/lib_kernels_maths.h~

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