src/snake2D_gpu.cu~

   1 /**
   2  * \file snake2D.c
   3  * \brief snake polygonale approche region sous hypothese Gaussienne ou Gamma.
   4  * \author NB - PhyTI
   5  * \version x.x
   6  * \date 20 decembre 2009
   7  *
   8  * traitement d'images en entiers 16 bits non signe : ppm
   9  * USAGE : SNAKE2D image.pgm 0 (ou 1)
  10  */
  11
  12 #include <stdio.h>
  13 #include <malloc.h>
  14 #include "structures.h"
  15 extern "C"{
  16 #include "lib_alloc.h"
  17 #include "lib_images.h"
  18 #include "lib_snake_common.h"
  19 #include "lib_math.h"
  20 #include "lib_gpu.h"
  21 #include "defines.h"
  22 #include "lib_snake_2_gpu.h"
  23 }
  24 #include "lib_kernels_maths.cu"
  25 #include "lib_kernels_contribs.cu"
  26
  27
  28
  29
  30 int main(int argc, char **argv)
  31 {
  32   /* declaration des variables */
  33   int ret ;
  34   int Prof ;        /* profondeur en octets */
  35   uint32 I_dim ;       /* hauteur de l'image */
  36   uint32 J_dim ;       /* largeur de l'image */
  37   int Nb_level ;    /* dynamique de l'image */
  38   char *File_name ;
  39
  40
  41   /* images */
  42   unsigned short **Image_in;
  43   struct timeval chrono, chrono2, chrono_all ;
  44
  45   /* lecture argument entree (basique!) */
  46   File_name = argv[3] ;
  47
  48   /* verif type image (pgm 8/16) */
  49   ret = type_image_ppm(&Prof, &I_dim, &J_dim, &Nb_level, File_name) ;
  50
  51   if ((ret == 0) | (Prof == 3))
  52     {
  53       printf("format non pris en charge ... exit\n") ;
  54       return(0) ;
  55     }
  56
  57   /* infos */
  58   printf("Image : %s\n", File_name) ;
  59   printf("lecture OK : %d\n", ret) ;
  60   printf("Image (%d x %d) pixels\n", I_dim, J_dim) ;
  61   printf("Dynamique : %d\n", Nb_level) ;
  62
  63   /* Allocation */
  64   Image_in = new_matrix_ushort(I_dim, J_dim) ;
  65
  66   /* chargement image d'entree */
  67   load_pgm2ushort(Image_in, I_dim, J_dim, Nb_level, File_name) ;
  68
  69   //POINTEURS VARIABLES MEMOIRE GLOBALE GPU
  70   unsigned short    * d_img ;   // image
  71   t_cumul_x * d_img_x ;         // images cumulees
  72   t_cumul_x2 * d_img_x2;        //
  73
  74   snake_node_gpu * d_snake ;    //image du snake CPU dans un tableau en gmem GPUe
  75
  76   int    * d_freemanDiDj ;            // table de correspondance [Di][Dj]->Freemans
  77   int    * d_codeNoeud ;              // table de correspondance [F_in][F_out]->codeNoeud
  78
  79   uint4  * d_positions ;              // positions de test autour des noeuds
  80
  81   uint2  * d_listes_pixels ;          // coordonnees des pixels des segments correspondants aux 8 posiionstest
  82   uint2  * d_liste_temp ;
  83   uint32 * d_nb_pix_max ;             // taille max des segments a tester
  84
  85   uint64 * d_contribs_segments_blocs ;// sommes des contribs pixels par blocs de calcul
  86   uint64 * d_contribs_segments ;      // contribs segments 1, x et x2
  87   uint64 * d_sompart ;                // vecteur de resultats intermediaires (sommes partielles = sommes par blocs)
  88
  89   int64  * d_stats, * d_stats_ref ;   // stats des positions de test, du snake sans les segments en test
  90   int64  * d_stats_snake;             // stats du snake + stats de l'image complete
  91   double * d_vrais, * d_vrais_snake ; // valeurs de la log-vraisemblance des positions de test, du snake
  92
  93   uint4  * d_freemans_centres ;       // valeurs des F_in, F_out et coord.
  94                                       // centres des 16 segments associes aux 8 positions de test
  95
  96   int    * d_codes_segments ;         // valeurs de codes des 16 segments
  97   bool   * d_move ;                   // nb de deplacement effectues lors d'une iteration
  98   int    * d_nb_nodes ;               // nb de noeuds du snake
  99
 100   snake_node_gpu * d_snake_tmp ;      // snake tampon pour l'etape d'ajout de noeuds
 101
 102   /* pointeurs sur mem CPU */
 103   int *h_nb_nodes = new int;          // image CPU du nb de noeud du snake
 104   snake_node_gpu  h_snake[MAX_NODES] ;
 105   uint2 h_listes_pixels[MAX_NODES*16*5] ;
 106   double h_vrais_snake, h_vrais_mem, h_vrais[8*MAX_NODES] ;              // image CPU de la log-vraisemblance
 107   bool * h_move = new bool[MAX_NODES];// image CPU du vecteur identifiant les noeuds qui ont bouge
 108   uint32 h_nb_pix_max, npixmax ;      // taille max des segments a tester : utile pour determiner les params d'execution
 109   int nnodes = 4 ;                    // 4 ou 40 pour l'instant
 110
 111
 112   /*allocation memoire GPU  */
 113   int retour_cmd = 0;
 114   cudaMalloc((void**) &d_nb_nodes, sizeof(int));
 115   cudaMalloc((void**) &d_sompart, MAX_NODES*256*16*sizeof(uint64));
 116   cudaMalloc((void**) &d_liste_temp, MAX_NODES*5*16*sizeof(uint2));
 117   retour_cmd = cudaMalloc((void**) &d_snake_tmp, MAX_NODES*sizeof(snake_node_gpu) );
 118   printf("RESULTAT ALLOC snake_tmp = %d\n", retour_cmd);
 119
 120
 121   /*init snake (positions/contribs/stats/freemans/centres/codes)*/
 122   cuda_init_img_cumul(Image_in, I_dim, J_dim, nnodes,
 123                                           &d_img, &d_img_x, &d_img_x2,
 124                                           &d_freemanDiDj, &d_codeNoeud,
 125                                           &d_snake, &d_nb_pix_max,
 126                                           &d_positions, &d_contribs_segments, &d_freemans_centres,
 127                                           &d_codes_segments, &d_stats_snake,
 128                                           &d_stats, &d_stats_ref, &d_vrais, &d_vrais_snake,
 129                                           &d_listes_pixels, &d_contribs_segments_blocs,
 130                                           &d_move
 131                                           );
 132
 133   /* debug : affichage snake */
 134   int Verbose = 1 ;
 135   int VERBOSE = 1 ;
 136   int Display = 1 ;
 137
 138   /* debug : tests fonctions*/
 139   int H= I_dim, L= J_dim ;
 140
 141   //snake_node * h_snake_ll;
 142   uint64 h_stats_snake[6];
 143   //gpu2snake(d_snake, &h_snake_ll, nnodes);
 144
 145
 146   // variables de debug
 147   int nb_move, nb_newnode, iter, i ;
 148   int nb_move_total=0, nb_test_total=0 ;
 149   int NB_iter_max = atoi(argv[1]);
 150   int dist = (I_dim+J_dim)/5;
 151   int Pas = atoi(argv[2]) ;                           // dist entre la position actuelle et les positions de test
 152   int Dist_min_entre_noeud = 4*Pas ;
 153   int bs, nblocs_seg, tpb, bps ;                  // nb de threads par blocs pour l'execution des kernels, nb de blocs de threads par segment a tester
 154   dim3 threads, grid ;                            // params d'execution des kernels
 155   int n_interval ;                                // nombre d'intervalles Na--Nx--Nb concernes
 156   int taille_smem ;                               // quantite de shared memory allouee pour le calcul des contribs des segments de test
 157   bool pairs = true ;                             // mouvement des noeuds pairs/impairs
 158
 159   if (Verbose) {
 160         printf("nb noeuds : %d\n", nnodes) ;
 161         tic(&chrono_all, NULL) ;
 162   }
 163
 164   for (iter=1; (iter<=NB_iter_max)&&(Pas>0); iter++, Pas>>=1)
 165     {
 166
 167       if (VERBOSE)
 168                 {
 169                   cudaMemcpy( &h_vrais_snake, d_vrais_snake, sizeof(double), cudaMemcpyDeviceToHost);
 170                   printf("\n#%d : pas %d pixels, LV = %lf \n", iter, Pas, h_vrais_snake) ;
 171                   tic(&chrono, NULL) ;
 172                 }
 173           // DEBUT MOVE SNAKE
 174           do {
 175
 176                 //memorisation precedente LV
 177                 h_vrais_mem = h_vrais_snake ;
 178                 // calcul stats sans les paires de segments a bouger
 179                 soustrait_aux_stats_2N_segments_noeud<<< nnodes , 1 >>>(d_snake, d_stats_snake, d_stats_ref,
 180                                                                         d_img_x, d_img_x2,
 181                                                                         d_codeNoeud, J_dim
 182                                                                         );
 183
 184                 // calcul des coordonnées de toutes les positions possibles des noeud a l'etape N+1
 185                 liste_positions_a_tester<<<nnodes, 8>>>(d_snake, d_positions, d_nb_pix_max, Pas, nnodes, I_dim, J_dim) ;
 186
 187                 // recupere la taille maxi des segments
 188                 cudaMemcpy( &h_nb_pix_max, d_nb_pix_max, sizeof(uint32), cudaMemcpyDeviceToHost) ;
 189
 190                 // determination des parametres des kernels
 191                 bs = nextPow2(h_nb_pix_max) ;
 192                 if (bs>=BSMAX) bs = BSMAX ; //  /!\ le kernel <<< calcul_contrib...>>> ne supporte pas un bs>256 a cause de la shared-mem nécessaire
 193                 if (bs<32) bs = 32 ;
 194                 nblocs_seg = (h_nb_pix_max+bs-1)/bs ;
 195
 196                 pairs = false ;
 197                 n_interval = nnodes/2 + pairs*(nnodes%2) ;
 198                 taille_smem =  CFI(bs)*sizeof(tcontribs) ;
 199                 threads = dim3(bs,1,1) ;
 200                 grid = dim3( n_interval*16*nblocs_seg ,1,1) ;
 201
 202                   //calcul listes pix + contrib partielles + freemans + centres
 203                   calcul_contribs_segments_blocs_full<<< grid , threads, taille_smem >>>( d_snake, nnodes, d_positions, h_nb_pix_max,
 204                                                                                           d_img_x, d_img_x2, d_codes_segments,
 205                                                                                           J_dim, d_listes_pixels, d_contribs_segments_blocs,
 206                                                                                                                                                                   pairs);
 207                   /*debug*/
 208                   cudaMemcpy( h_listes_pixels, d_listes_pixels, 16*5*n_interval*sizeof(uint2), cudaMemcpyDeviceToHost);
 209                   for(int inter=0; inter < 3; inter++){
 210                         for(int seg=0; seg<16; seg++){
 211                           printf(" intervalle %d segment %d : (%d,%d)-(%d,%d)-(%d,%d)-(%d,%d)-(%d,%d)\n",
 212                                          inter, seg, h_listes_pixels[5*(16*inter+seg)].x,h_listes_pixels[5*(16*inter+seg)].y,h_listes_pixels[5*(16*inter+seg)+1].x,
 213                                          h_listes_pixels[5*(16*inter+seg)+1].y,h_listes_pixels[5*(16*inter+seg)+2].x,h_listes_pixels[5*(16*inter+seg)+2].y,
 214                                 h_listes_pixels[5*(16*inter+seg)+3].x,h_listes_pixels[5*(16*inter+seg)+3].y,h_listes_pixels[5*(16*inter+seg)+4].x,
 215                                 h_listes_pixels[5*(16*inter+seg)+4].y);
 216                         }
 217                   }
 218                   /*fin*/
 219           calcul_freemans_centre<<<n_interval, 16>>>( d_listes_pixels,  d_freemanDiDj, d_freemans_centres);
 220                   //printf("EXEC impairs : %d max pix - %d intervalles => %d blocs de %d threads - %d octets de smem\n", h_nb_pix_max, n_interval, grid.x, threads.x, taille_smem);
 221                   //sommes des contribs partielles -> contribs segments
 222                   somsom_full<<< 16*n_interval , 1>>>(d_contribs_segments_blocs, nnodes, nblocs_seg, d_contribs_segments) ;
 223
 224                   //calcul des stats associees a chaque position de test
 225                   calcul_stats_full<<< n_interval, 8 >>>(d_snake, nnodes, pairs, d_stats_snake, d_stats_ref, d_stats, d_contribs_segments,
 226                                                                                                  d_positions, d_codes_segments,  d_freemans_centres, d_codeNoeud,
 227                                                                                                  d_img_x, d_img_x2, I_dim, J_dim, d_vrais, d_vrais_snake, d_move);
 228
 229                   pairs = true ;
 230                   n_interval = nnodes/2 + pairs*(nnodes%2) ;
 231                   grid = dim3( n_interval*16*nblocs_seg ,1,1) ;
 232
 233                   //calcul listes pix + contrib partielles + freemans + centres
 234                   calcul_contribs_segments_blocs_full<<< grid , threads, taille_smem >>>( d_snake, nnodes, d_positions, h_nb_pix_max,
 235                                                                                           d_img_x, d_img_x2, d_codes_segments,
 236                                                                                           J_dim, d_listes_pixels, d_contribs_segments_blocs,
 237                                                                                           pairs);
 238                 calcul_freemans_centre<<<n_interval, 16>>>( d_listes_pixels, d_freemanDiDj, d_freemans_centres);
 239                 //printf("EXEC pairs : %d max pix - %d intervalles => %d blocs de %d threads - %d octets de smem\n", h_nb_pix_max, n_interval, grid.x, threads.x, taille_smem);
 240                 //sommes des contribs partielles -> contribs segments
 241                 somsom_full<<< 16*n_interval , 1>>>(d_contribs_segments_blocs, nnodes, nblocs_seg, d_contribs_segments) ;
 242
 243                 //calcul des stats associees a chaque position de test
 244                 calcul_stats_full<<< n_interval, 8 >>>(d_snake, nnodes, pairs, d_stats_snake, d_stats_ref, d_stats, d_contribs_segments,
 245                                                                                                  d_positions, d_codes_segments,  d_freemans_centres, d_codeNoeud,
 246                                                                                                  d_img_x, d_img_x2, I_dim, J_dim, d_vrais, d_vrais_snake, d_move);
 247
 248
 249                 //il faut recalculer les stats du snake apres modif
 250                 recalcul_stats_snake<<< 1 , 1  >>>(d_snake, nnodes, d_stats_snake, d_vrais_snake,
 251                                                                                    d_img_x, d_img_x2,
 252                                                                                    d_codeNoeud, J_dim
 253                                                                                    );
 254
 255                 cudaMemcpy( &h_vrais_snake, d_vrais_snake, sizeof(double), cudaMemcpyDeviceToHost);
 256                 //printf("iter %d apres recalcul du move LV = %lf - ",  iter, h_vrais_snake) ;
 257
 258                 nb_move = 0;
 259                 //recup move
 260                 cudaMemcpy( h_move, d_move, nnodes*sizeof(bool), cudaMemcpyDeviceToHost);
 261                 i = 0;
 262                 while (i<nnodes)
 263                   {
 264                         nb_move += (int)h_move[i];
 265                         i++;
 266                   }
 267
 268                 nb_move_total += nb_move ;
 269                 nb_test_total+= nnodes*8 ;
 270           } while ( nb_move && (h_vrais_snake < h_vrais_mem));
 271
 272           if ( iter < NB_iter_max ){
 273             // ajout de noeud et recalcul stats
 274             ajoute_noeuds<<< 1 , 1 >>>(d_snake, d_snake_tmp, nnodes, Dist_min_entre_noeud, d_nb_nodes );
 275             cudaMemcpy( h_nb_nodes, d_nb_nodes, sizeof(int), cudaMemcpyDeviceToHost);
 276             nnodes += (*h_nb_nodes) ;
 277
 278             cudaMemcpy( d_snake, d_snake_tmp, nnodes*sizeof(snake_node_gpu), cudaMemcpyDeviceToDevice);
 279
 280             npixmax = h_nb_pix_max ;
 281             tpb = nextPow2(npixmax) ;
 282             if (tpb >= BSMAX) tpb = BSMAX ;//  /!\ le kernel <<< calcul_contrib...>>> ne supporte pas un bs>256 a cause de la shared-mem nécessaire
 283             if (tpb < 32 ) tpb = 32 ;
 284             bps = (npixmax+tpb-1)/tpb ;
 285
 286             recalcul_contribs_segments_snake<<< nnodes*bps, tpb, CFI(tpb)*sizeof(tcontribs)>>>(d_snake, nnodes,
 287                                                                                                 d_img_x, d_img_x2,
 288                                                                                                 J_dim, d_liste_temp, d_sompart );
 289
 290             recalcul_freemans_centre<<<nnodes, 1>>>(d_snake, d_liste_temp, d_freemanDiDj);
 291             resomsom_snake<<< nnodes , 1 >>>(d_sompart, nnodes, bps, d_snake);
 292
 293             recalcul_stats_snake<<< 1 , 1 >>>(d_snake, nnodes, d_stats_snake, d_vrais_snake,
 294                                                   d_img_x, d_img_x2,
 295                                                   d_codeNoeud, J_dim
 296                                                   );
 297
 298                 if (*h_nb_nodes == 0) break ;
 299                 cudaMemcpy( &h_vrais_snake, d_vrais_snake, sizeof(double), cudaMemcpyDeviceToHost);
 300                 printf("\niter %d LV apres ajout noeuds = %lf \n ", iter, h_vrais_snake) ;
 301           }
 302
 303       if (VERBOSE)
 304         {
 305           toc(chrono, "temps sequence move");
 306
 307           printf("nb deplacements    : %d\n", nb_move) ;
 308           printf("nb deplacements total/test   : %d/%d\n", nb_move_total, nb_test_total) ;
 309           printf("nb nouveaux noeuds : %d (total: %d)\n", *h_nb_nodes, nnodes) ;
 310           printf("\nlongueur de codage de gl : %lf  \n", h_vrais_snake) ;
 311         }
 312     }
 313
 314   if (Verbose) {
 315         toc(chrono_all, "temps move mv") ;
 316         cudaMemcpy( h_stats_snake, d_stats_snake, 6*sizeof(uint64), cudaMemcpyDeviceToHost);
 317         cudaMemcpy( &h_vrais_snake, d_vrais_snake, sizeof(double), cudaMemcpyDeviceToHost);
 318         printf("\nFIN : longueur de codage de gl : %lf  (%d)\n", h_vrais_snake, h_stats_snake[0]) ;
 319         printf("nb noeuds : %d, nb_iter : %d\n", nnodes, iter-1) ;
 320         printf("nb deplacements total/test   : %d/%d\n", nb_move_total, nb_test_total) ;
 321   }
 322
 323
 324   if (Display) {
 325         /* old fashion way to draw the snake
 326         gpu2snake(d_snake, &h_snake_ll, nnodes);
 327         uint32 * Liste_pixel_segment = new uint32[I_dim+J_dim];
 328         affiche_snake_ushort(Image_in, h_snake_ll, 255, 0, Liste_pixel_segment) ;
 329         delete Liste_pixel_segment ;
 330         delete h_snake_ll;
 331         */
 332         cudaMemcpy( h_snake, d_snake, nnodes*sizeof(snake_node_gpu), cudaMemcpyDeviceToHost);
 333         //affiche coordonnees
 334         for (int node=0; node<nnodes; node++){
 335             printf("NODE %d  %d  %d \n", node, h_snake[node].posi, h_snake[node].posj);
 336         }
 337         // draw only the nodes with + marks
 338         //dessine_snake(h_snake, nnodes, Image_in, 10);
 339         //imagesc_ushort(Image_in, I_dim, J_dim) ;
 340   }
 341   cudaFree(d_img);
 342   cudaFree(d_img_x);
 343   cudaFree(d_img_x2);
 344   cudaFree(d_freemanDiDj);
 345   cudaFree(d_codeNoeud);
 346   cudaFree(d_snake);
 347   cudaFree(d_nb_pix_max);
 348   cudaFree(d_positions);
 349   cudaFree(d_contribs_segments);
 350   cudaFree(d_freemans_centres);
 351   cudaFree(d_codes_segments);
 352   cudaFree(d_stats_snake);
 353   cudaFree(d_stats);
 354   cudaFree(d_stats_ref);
 355   cudaFree(d_vrais);
 356   cudaFree(d_vrais_snake);
 357   cudaFree(d_listes_pixels);
 358   cudaFree(d_contribs_segments_blocs);
 359   cudaFree(d_move);
 360   return 1 ;
 361 }
 362