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, 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   double h_vrais_snake, h_vrais_mem ; // image CPU de la log-vraisemblance
 106   bool * h_move = new bool[MAX_NODES];// image CPU du vecteur identifiant les noeuds qui ont bouge
 107   uint32 h_nb_pix_max, npixmax ;      // taille max des segments a tester : utile pour determiner les params d'execution
 108   int nnodes = 4 ;
 109
 110   /*allocation memoire GPU  */
 111   cudaMalloc((void**) &d_nb_nodes, sizeof(int));
 112   cudaMalloc((void**) &d_sompart, MAX_NODES*256*16*sizeof(uint64));
 113   cudaMalloc((void**) &d_liste_temp, MAX_NODES*5*16*sizeof(uint2));
 114   cudaMalloc((void**) &d_snake_tmp, MAX_NODES*sizeof(snake_node_gpu) );
 115
 116   /*init snake (positions/contribs/stats/freemans/centres/codes)*/
 117   cuda_init_img_cumul(Image_in, I_dim, J_dim, nnodes,
 118                                           &d_img, &d_img_x, &d_img_x2,
 119                                           &d_freemanDiDj, &d_codeNoeud,
 120                                           &d_snake, &d_nb_pix_max,
 121                                           &d_positions, &d_contribs_segments, &d_freemans_centres,
 122                                           &d_codes_segments, &d_stats_snake,
 123                                           &d_stats, &d_stats_ref, &d_vrais, &d_vrais_snake,
 124                                           &d_listes_pixels, &d_contribs_segments_blocs,
 125                                           &d_move
 126                                           );
 127
 128   /* debug : affichage snake */
 129   int Verbose = 1 ;
 130   int VERBOSE = 1 ;
 131   int Display = 1 ;
 132
 133   uint64 h_stats_snake[6];
 134
 135   // variables de debug
 136   int nb_move, iter, i ;
 137   int nb_move_total=0, nb_test_total=0 ;
 138   int NB_iter_max = atoi(argv[1]);
 139   int Pas = atoi(argv[2]) ;                       // distance entre la position actuelle et les positions de test
 140   int Dist_min_entre_noeud = 4*Pas ;
 141   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
 142   dim3 threads, grid ;                            // params d'execution des kernels
 143   int n_interval ;                                // nombre d'intervalles Na--Nx--Nb concernes
 144   int taille_smem ;                               // quantite de shared memory allouee pour le calcul des contribs des segments de test
 145   bool pairs = true ;                             // mouvement des noeuds pairs/impairs
 146
 147   if (Verbose) {
 148         printf("nb noeuds : %d\n", nnodes) ;
 149         tic(&chrono_all, NULL) ;
 150   }
 151
 152   for (iter=1; (iter<=NB_iter_max)&&(Pas>0); iter++, Pas>>=1)
 153     {
 154       if (VERBOSE)
 155                 {
 156                   cudaMemcpy( &h_vrais_snake, d_vrais_snake, sizeof(double), cudaMemcpyDeviceToHost);
 157                   printf("\n#%d : pas %d pixels, LV = %lf \n", iter, Pas, h_vrais_snake) ;
 158                   tic(&chrono, NULL) ;
 159                 }
 160           // DEBUT MOVE SNAKE
 161           do {
 162
 163                 //memorisation precedente LV
 164                 h_vrais_mem = h_vrais_snake ;
 165                 // calcul stats sans les paires de segments a bouger
 166                 soustrait_aux_stats_2N_segments_noeud<<< nnodes , 1 >>>(d_snake, d_stats_snake, d_stats_ref,
 167                                                                         d_img_x, d_img_x2,
 168                                                                         d_codeNoeud, J_dim
 169                                                                         );
 170
 171                 // calcul des coordonnées de toutes les positions possibles des noeud a l'etape N+1
 172                 liste_positions_a_tester<<<nnodes, 8>>>(d_snake, d_positions, d_nb_pix_max, Pas, nnodes, I_dim, J_dim) ;
 173
 174                 // recupere la taille maxi des segments
 175                 cudaMemcpy( &h_nb_pix_max, d_nb_pix_max, sizeof(uint32), cudaMemcpyDeviceToHost) ;
 176
 177                 // determination des parametres des kernels
 178                 bs = nextPow2(h_nb_pix_max) ;
 179                 if (bs>=BSMAX) bs = BSMAX ; //  /!\ le kernel <<< calcul_contrib...>>> ne supporte pas un bs>256 a cause de la shared-mem nécessaire
 180                 if (bs<32) bs = 32 ;
 181                 nblocs_seg = (h_nb_pix_max+bs-1)/bs ;
 182
 183                 pairs = false ;
 184                 n_interval = nnodes/2 + pairs*(nnodes%2) ;
 185                 taille_smem =  CFI(bs)*sizeof(tcontribs) ;
 186                 threads = dim3(bs,1,1) ;
 187                 grid = dim3( n_interval*16*nblocs_seg ,1,1) ;
 188
 189                   //calcul listes pix + contrib partielles + freemans + centres
 190                   calcul_contribs_segments_blocs_full<<< grid , threads, taille_smem >>>( d_snake, nnodes, d_positions, h_nb_pix_max,
 191                                                                                           d_img_x, d_img_x2, d_codes_segments,
 192                                                                                           J_dim, d_listes_pixels, d_contribs_segments_blocs,
 193                                                                                                                                                                   pairs);
 194
 195           calcul_freemans_centre<<<n_interval, 16>>>( d_listes_pixels,  d_freemanDiDj, d_freemans_centres);
 196                   //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);
 197                   //sommes des contribs partielles -> contribs segments
 198                   somsom_full<<< 16*n_interval , 1>>>(d_contribs_segments_blocs, nnodes, nblocs_seg, d_contribs_segments) ;
 199
 200                   //calcul des stats associees a chaque position de test
 201                   calcul_stats_full<<< n_interval, 8 >>>(d_snake, nnodes, pairs, d_stats_snake, d_stats_ref, d_stats, d_contribs_segments,
 202                                                                                                  d_positions, d_codes_segments,  d_freemans_centres, d_codeNoeud,
 203                                                                                                  d_img_x, d_img_x2, I_dim, J_dim, d_vrais, d_vrais_snake, d_move);
 204
 205                   pairs = true ;
 206                   n_interval = nnodes/2 + pairs*(nnodes%2) ;
 207                   grid = dim3( n_interval*16*nblocs_seg ,1,1) ;
 208
 209                   //calcul listes pix + contrib partielles + freemans + centres
 210                   calcul_contribs_segments_blocs_full<<< grid , threads, taille_smem >>>( d_snake, nnodes, d_positions, h_nb_pix_max,
 211                                                                                           d_img_x, d_img_x2, d_codes_segments,
 212                                                                                           J_dim, d_listes_pixels, d_contribs_segments_blocs,
 213                                                                                           pairs);
 214                   calcul_freemans_centre<<<n_interval, 16>>>( d_listes_pixels, d_freemanDiDj, d_freemans_centres);
 215                   //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);
 216                 //sommes des contribs partielles -> contribs segments
 217                   somsom_full<<< 16*n_interval , 1>>>(d_contribs_segments_blocs, nnodes, nblocs_seg, d_contribs_segments) ;
 218
 219                   //calcul des stats associees a chaque position de test
 220                   calcul_stats_full<<< n_interval, 8 >>>(d_snake, nnodes, pairs, d_stats_snake, d_stats_ref, d_stats, d_contribs_segments,
 221                                                                                                  d_positions, d_codes_segments,  d_freemans_centres, d_codeNoeud,
 222                                                                                                  d_img_x, d_img_x2, I_dim, J_dim, d_vrais, d_vrais_snake, d_move);
 223
 224
 225                 //il faut recalculer les stats du snake apres modif
 226                 recalcul_stats_snake<<< 1 , 1  >>>(d_snake, nnodes, d_stats_snake, d_vrais_snake,
 227                                                                                    d_img_x, d_img_x2,
 228                                                                                    d_codeNoeud, J_dim
 229                                                                                    );
 230
 231                 cudaMemcpy( &h_vrais_snake, d_vrais_snake, sizeof(double), cudaMemcpyDeviceToHost);
 232                 //printf("iter %d apres recalcul du move LV = %lf - ",  iter, h_vrais_snake) ;
 233
 234                 nb_move = 0;
 235                 //recup move
 236                 cudaMemcpy( h_move, d_move, nnodes*sizeof(bool), cudaMemcpyDeviceToHost);
 237                 i = 0;
 238                 while (i<nnodes)
 239                   {
 240                         nb_move += (int)h_move[i];
 241                         i++;
 242                   }
 243
 244                 nb_move_total += nb_move ;
 245                 nb_test_total+= nnodes*8 ;
 246           } while ( nb_move && (h_vrais_snake < h_vrais_mem));
 247
 248           if ( iter < NB_iter_max ){
 249             // ajout de noeuds
 250             ajoute_noeuds<<< 1 , 1 >>>(d_snake, d_snake_tmp, nnodes, Dist_min_entre_noeud, d_nb_nodes );
 251                 //recup nb de nouveaux noeuds
 252                 cudaMemcpy( h_nb_nodes, d_nb_nodes, sizeof(int), cudaMemcpyDeviceToHost);
 253                 //mise a jour nb de noeuds
 254             nnodes += (*h_nb_nodes) ;
 255
 256                 //parametres d'execution des kernels pour le recalcul des contribs et stats du snake
 257             npixmax = h_nb_pix_max ;
 258             tpb = nextPow2(npixmax) ;
 259             if (tpb >= BSMAX) tpb = BSMAX ;//  /!\ le kernel <<< calcul_contrib...>>> ne supporte pas un bs>BSMAX a cause de la shared-mem nécessaire
 260             if (tpb < 32 ) tpb = 32 ;
 261             bps = (npixmax+tpb-1)/tpb ;
 262
 263                 //calcul sommes partielles des contribs + codes segments
 264             recalcul_contribs_segments_snake<<< nnodes*bps, tpb, CFI(tpb)*sizeof(tcontribs)>>>(d_snake, nnodes,
 265                                                                                                                                                                                    d_img_x, d_img_x2,
 266                                                                                                                                                                                    J_dim, d_liste_temp, d_sompart );
 267                 //calcul des freemans et des centres a partir des 5 points stockes par segment dans 'd_liste_temp'
 268             recalcul_freemans_centre<<<nnodes, 1>>>(d_snake, d_liste_temp, d_freemanDiDj);
 269                 //somme des sommes partielles
 270             resomsom_snake<<< nnodes , 1 >>>(d_sompart, nnodes, bps, d_snake);
 271                 //calcul des stats
 272             recalcul_stats_snake<<< 1 , 1 >>>(d_snake, nnodes, d_stats_snake, d_vrais_snake,
 273                                                                                   d_img_x, d_img_x2,
 274                                                                                   d_codeNoeud, J_dim
 275                                                                                   );
 276             //tant que l'on peut ajouter des noeuds
 277                 if (*h_nb_nodes == 0) break ;
 278                 //recup LogVraisemblance
 279                 cudaMemcpy( &h_vrais_snake, d_vrais_snake, sizeof(double), cudaMemcpyDeviceToHost);
 280
 281           }
 282
 283       if (VERBOSE)
 284         {
 285           toc(chrono, "temps sequence move");
 286           printf("nb deplacements    : %d\n", nb_move) ;
 287           printf("nb deplacements total/test   : %d/%d\n", nb_move_total, nb_test_total) ;
 288           printf("nb nouveaux noeuds : %d (total: %d)\n", *h_nb_nodes, nnodes) ;
 289           printf("\nlongueur de codage de gl : %lf  \n", h_vrais_snake) ;
 290         }
 291     }
 292
 293   if (Verbose) {
 294         toc(chrono_all, "temps move mv") ;
 295         cudaMemcpy( h_stats_snake, d_stats_snake, 6*sizeof(uint64), cudaMemcpyDeviceToHost);
 296         cudaMemcpy( &h_vrais_snake, d_vrais_snake, sizeof(double), cudaMemcpyDeviceToHost);
 297         printf("\nFIN : longueur de codage gl : %lf  (%d)\n", h_vrais_snake, h_stats_snake[0]) ;
 298         printf("nb noeuds : %d, nb_iter : %d\n", nnodes, iter-1) ;
 299         printf("nb deplacements total/test   : %d/%d\n", nb_move_total, nb_test_total) ;
 300   }
 301
 302
 303   if (Display) {
 304         cudaMemcpy( h_snake, d_snake, nnodes*sizeof(snake_node_gpu), cudaMemcpyDeviceToHost);
 305         //affiche coordonnees
 306         for (int node=0; node<nnodes; node++){
 307             printf("NODE %d =  ( %d , %d ) \n", node, h_snake[node].posi, h_snake[node].posj);
 308         }
 309         // dessine seulement les  nodes avec des +
 310         dessine_snake(h_snake, nnodes, Image_in, 10);
 311         imagesc_ushort(Image_in, I_dim, J_dim) ;
 312   }
 313   cudaFree(d_img);
 314   cudaFree(d_img_x);
 315   cudaFree(d_img_x2);
 316   cudaFree(d_freemanDiDj);
 317   cudaFree(d_codeNoeud);
 318   cudaFree(d_snake);
 319   cudaFree(d_nb_pix_max);
 320   cudaFree(d_positions);
 321   cudaFree(d_contribs_segments);
 322   cudaFree(d_freemans_centres);
 323   cudaFree(d_codes_segments);
 324   cudaFree(d_stats_snake);
 325   cudaFree(d_stats);
 326   cudaFree(d_stats_ref);
 327   cudaFree(d_vrais);
 328   cudaFree(d_vrais_snake);
 329   cudaFree(d_listes_pixels);
 330   cudaFree(d_contribs_segments_blocs);
 331   cudaFree(d_move);
 332   return 0 ;
 333 }
 334