]> AND Private Git Repository - snake_gpu.git/blob - src/snake2D_gpu.cu~
Logo AND Algorithmique Numérique Distribuée

Private GIT Repository
Fin de test multisnake sur iter 1
[snake_gpu.git] / 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