+#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];
+ }
+}
+