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