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