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