Fin de test multisnake sur iter 1

[snake_gpu.git] / src / lib_kernel_snake_2_gpu.cu

diff --git a/src/lib_kernel_snake_2_gpu.cu b/src/lib_kernel_snake_2_gpu.cu
index b807cb77e519efa15cfbdaf9a6c80a2367efb455..3796e6569c7e152a73e4bc1f96885610a07f06a9 100644 (file)
--- a/src/lib_kernel_snake_2_gpu.cu
+++ b/src/lib_kernel_snake_2_gpu.cu
@@ -33,7 +33,7 @@ __global__ void genere_snake_rectangle_4nodes_gpu(snake_node_gpu * d_snake, int
   }
 }
 
-__global__ void genere_diagos_rectangle(uint4 * d_diagos, int h, int l, int q){
+__global__ void genere_diagos_rectangle(uint4 * d_diagos, int h, int l, int q, int * n_diagos){
   int inci = h/q;
   int incj = l/q;
   int iM,jM, iN, jN ;
@@ -53,6 +53,7 @@ __global__ void genere_diagos_rectangle(uint4 * d_diagos, int h, int l, int q){
                }
          }
        }
+       *n_diagos = --idxDiago ;
 }
 
 __global__ void genere_snake_rectangle_Nnodes_gpu(snake_node_gpu * d_snake, int dist_bords, int i_dim, int j_dim){
@@ -466,3 +467,88 @@ __global__ void calcul_stats_snake(snake_node_gpu * d_snake, int  nnodes, int64
   *vrais_min = codage_gl_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]);
 }
+
+// kernel d'init rectangulaire au niveau snake
+// un block par point de base et un point opposé par thread du bloc
+
+__global__ void calcul_contribs_snake4(snake_node_gpu * snake, t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2, int h, int l,
+                                                                          int64 * sums, t_rectangle_snake * critere)
+{
+  // nb de diagonales testees par bloc (ie. par point de base NO)
+  int blockSize = blockDim.x ;
+  // indice du second point de chaque diagonale (=Opposite Point, = point SE)
+  int OPib = threadIdx.x ;              
+  // coordonnees de chaque point de base (NO)
+  int BPi = blockIdx.x ;
+  int BPj = blockIdx.y ;
+  //coordonnees de chaque Opposite Point (SE)
+  double incThread = ((h-BPi)*(l-BPj) + blockSize -1)/(double)blockSize ;
+  int OPi = (double)(OPib*incThread)/(l - BPj) ;
+  int OPj = OPib*incThread - OPi*(l-BPj) ;
+  OPi += BPi ;
+  OPj += BPj ;
+  if (OPi >= h) OPi = h-1 ;
+  if (OPj >= l) OPj = l-1 ;
+  //indices des pixels dans les images cumulees
+  int posG, posD;
+  //contrib 1 du snake
+  int C1 = (OPi - BPi)*(OPj - BPj) ; 
+ 
+  //pour stocker contribs de chaque snake d'un block
+  //TODO on peut utiliser une structure restreinte (sans le c1)
+  extern __shared__ tcontribs scumuls[]; 
+  
+  //calcul contribs du snake
+  scumuls[CFI(OPib)].cx = 0;
+  scumuls[CFI(OPib)].cx2 = 0;
+  for (int k=BPi ; k < OPi ; k++)
+       {
+         posG = (BPi+k)*l + BPj ;
+         posD = posG - BPj + OPj ;
+         scumuls[CFI(OPib)].cx  += (cumul_x[ posD ] - cumul_x[ posG ]) ;
+         scumuls[CFI(OPib)].cx2 += (cumul_x2[ posD ] - cumul_x2[ posG ]);
+  } 
+  
+  //calcul de critère pour chaque snake
+  uint64 stat_sum_xe ;  /* somme des xn region exterieure */
+  uint32 ne ;           /* nombre de pixel region exterieure */
+  double sigi2, sige2;  /* variance region interieure et exterieure */ 
+  int index_crit;
+  
+  /* variance des valeurs des niveaux de gris a l'interieur du snake */
+  sigi2 = 
+    ((double)scumuls[CFI(OPib)].cx2/(double)C1) - 
+    ((double)scumuls[CFI(OPib)].cx/(uint64)C1)*((double)scumuls[CFI(OPib)].cx/(uint64)C1) ;
+
+  /* variance des valeurs des niveaux de gris a l'exterieur du snake */
+  
+  ne = sums[3] - C1 ;
+  stat_sum_xe = sums[4] - scumuls[CFI(OPib)].cx ;
+  sige2 =
+    ((double)sums[5]-scumuls[CFI(OPib)].cx2)/(double)ne - 
+    ((double)stat_sum_xe/(uint64)ne)*((double)stat_sum_xe/(uint64)ne) ;
+  
+  index_crit = blockSize*(BPi*gridDim.y + BPj) + OPib ;
+
+  
+  if ((sigi2 > 0)|(sige2 > 0))
+       {
+         critere[ index_crit ].bpi = BPi;
+         critere[ index_crit ].bpj = BPj;
+         critere[ index_crit ].opi = OPi;
+         critere[ index_crit ].opj = OPj;
+         critere[ index_crit ].crit = 0.5*((double)C1*log(sigi2) + (double)ne*log(sige2)) ;
+       }
+  else
+       {
+         critere[ index_crit ].bpi = BPi;
+         critere[ index_crit ].bpj = BPj;
+         critere[ index_crit ].opi = OPi;
+         critere[ index_crit ].opj = OPj;
+         critere[ index_crit ].crit = -1 ;
+       }
+  
+  // identification meilleur snake du bloc
+  // laissé au CPU pour test mais le principe de ce kernel n'est pas efficace.
+    
+}