]> AND Private Git Repository - lniv_gpu.git/blobdiff - levelines_kernels.cu
Logo AND Algorithmique Numérique Distribuée

Private GIT Repository
version operationnelle avec chemins parametrables
[lniv_gpu.git] / levelines_kernels.cu
index 491c600fbbea0f7aae989d02ba208447ffc216c5..6838e8322a7e2e786910c26cbaa0a5132f5ccd09 100644 (file)
-
-// chemins des lignes de niveaux
-// longueur = 4 pixels
-// une ligne = un chemin
-
+/************************************************************************
+ * chemins des lignes de niveaux pour la version à chemins constants
+ * Ne sont conservés que pour comparaison GPU/CPU -> à faire disparaître
+ * longueur = 4 pixels
+ * une ligne = un chemin
+ ************************************************************************/
 __constant__ int pathDi[PSIZE_I][PSIZE_J-1] =
   {
     // Q1
 __constant__ int pathDi[PSIZE_I][PSIZE_J-1] =
   {
     // Q1
-       {  -1, -1, -1},       // 
-       {  -1, -1, -1},       //  
-       {  -1, -1, -1},       // 
-       {  -1, -1, -1},       // 
-       {  -1,  0,  1},       // 
-       {   0, -1,  0},
+       {  -1, -1, -1},       // 90
+       {  -1, -1, -1},       // 75
+       {  -1, -1, -1},       // 60
+       {  -1, -1, -1},       // 45
+       {  -1,  0, -1},       // 30
+       {   0, -1,  0},       // 15
        // Q4
        // Q4
-       {   0,  0,  0},       // 
-       {   0,  1,  1},       //  
-       {   1,  0,  1},       // 
-       {   1,  1,  1},       // 
-       {   1,  1,  1},       // 
-       {   1,  1,  1},
+       {   0,  0,  0},       // 0
+       {   0,  1,  0},       // 345 
+       {   1,  0,  1},       // 330
+       {   1,  1,  1},       // 315
+       {   1,  1,  1},       // 300
+       {   1,  1,  1},       // 285
        // Q3
        // Q3
-       {   1,  1,  1},       // 
-       {   1,  1,  1},       //  
-       {   1,  1,  1},       // 
-       {   1,  1,  1},       // 
-       {   1,  0, -1},       // 
-       {   0,  1,  0},
+       {   1,  1,  1},       // 270
+       {   1,  1,  1},       // 255 
+       {   1,  1,  1},       // 240
+       {   1,  1,  1},       // 225
+       {   1,  0,  1},       // 210
+       {   0,  1,  0},       // 195
        // Q2
        // Q2
-       {   0,  0,  0},       // 
-       {   0, -1,  0},       //  
-       {  -1,  0, -1},       // 
-       {  -1, -1, -1},       // 
-       {  -1, -1, -1},       // 
-       {  -1, -1, -1}
+       {   0,  0,  0},       // 180
+       {   0, -1,  0},       // 165 
+       {  -1,  0, -1},       // 150
+       {  -1, -1, -1},       // 135
+       {  -1, -1, -1},       // 120
+       {  -1, -1, -1}        // 105
   } ;     // 
   
 __constant__ int pathDj[PSIZE_I][PSIZE_J-1] =
   {
        // Q1
   } ;     // 
   
 __constant__ int pathDj[PSIZE_I][PSIZE_J-1] =
   {
        // Q1
-       {  0,  0,  0},       // 
-       {  0,  1,  0},
-       {  1,  0,  1},  
-       {  1,  1,  1},  
-       {  1,  1,  1},  
-       {  1,  1,  1},
+       {  0,  0,  0},       // 90
+       {  0,  1,  0},       // 75
+       {  1,  0,  1},       // 60
+       {  1,  1,  1},       // 45
+       {  1,  1,  1},       // 30 
+       {  1,  1,  1},       // 15
        // Q4
        // Q4
-       {  1,  1,  1},       // 
-       {  1,  1,  1},
-       {  1,  1,  1},  
-       {  1,  1,  1},  
-       {  1,  0, -1},  
-       {  0,  1,  0},
-       // Q3
-       {  0,  0,  0},       // 
-       {  0, -1,  0},
-       { -1,  0, -1},  
-       { -1, -1, -1},  
-       { -1, -1, -1},  
-       { -1, -1, -1},
+       {  1,  1,  1},       // 0
+       {  1,  1,  1},       // 345
+       {  1,  1,  1},       // 330
+       {  1,  1,  1},       // 315  
+       {  1,  0,  1},       // 300
+       {  0,  1,  0},       // 285
+       // Q3 
+       {  0,  0,  0},       // 270
+       {  0, -1,  0},       // 255
+       { -1,  0, -1},       // 240
+       { -1, -1, -1},       // 225
+       { -1, -1, -1},       // 210
+       { -1, -1, -1},       // 195
        // Q2
        // Q2
-       { -1, -1, -1},       // 
-       { -1, -1, -1},
-       { -1, -1, -1},  
-       { -1, -1, -1},  
-       { -1,  0,  1},  
-       {  0, -1,  0}
+       { -1, -1, -1},       // 180
+       { -1, -1, -1},       // 165 
+       { -1, -1, -1},       // 150
+       { -1, -1, -1},       // 135
+       { -1,  0, -1},       // 120 
+       {  0, -1,  0}        // 105
   } ;     
 
   } ;     
 
+
+
+// valeurs des tangentes des angles de base pour la génération initiale des chemins
+// pour la version à chemins de longueur paramétrable
+__constant__ float tangente[] = {0.000, 0.268, 0.577, 1.000} ;
   
   
-// declare texture reference for 2D int texture
+// declarations des textures 
 texture<int, 2, cudaReadModeElementType> tex_img_in ;
 texture<int, 2, cudaReadModeElementType> tex_img_estim ;
 texture<int, 2, cudaReadModeElementType> tex_img_lniv ;
 texture<int, 2, cudaReadModeElementType> tex_img_in ;
 texture<int, 2, cudaReadModeElementType> tex_img_estim ;
 texture<int, 2, cudaReadModeElementType> tex_img_lniv ;
+texture<int2, 2, cudaReadModeElementType> tex_paths ;
+
+
+
+/**
+ *
+ * \brief calcule les chemins
+ * \author NB - PhyTI, modifié by zulu pour adaptater aux chemins paramétrables
+ *
+ * \param[in] r         longueur des chemins
+ *
+ * \param[out] d_paths  matrice des déplacements relatifs (chemins)
+ *
+ * Cette fonction utilise le tableau constant des tangentes des angles
+ * considérés pour le calcul de chemins (float tangente[]).
+ * 
+ */
+__global__ void kernel_calcul_paths( int2 * d_paths, unsigned int r){
+
+  unsigned int idpath = 0 ;
+  int ic, jc, iprec, jprec ;
+  float offset = 0.5 ;
+  unsigned int basepath = 0 ;
+
+  // Q1 inf
+  for (int a=0 ; a< 4 ; a++){        // les 4 angles 0,15,30 et 45
+       for (int p=0 ; p< r ; p++){      // les r points
+         ic = r-1 - floor(tangente[a]*p + offset) ;
+         if ( p > 0 ){
+               d_paths[idpath*(r-1)+p-1].x = ic - iprec ;
+               d_paths[idpath*(r-1)+p-1].y = 1 ;
+         }
+         iprec = ic ;
+       }
+       idpath++ ;
+  }
+  // Q1 sup
+  for (int a=2 ; a>0 ; a--){         // les 2 angles 60 et 75 
+       for (int p=0 ; p< r ; p++){      // les r points
+         jc = floor(tangente[a]*p + offset) ; 
+         if ( p > 0 ){
+               d_paths[idpath*(r-1)+p-1].x = -1 ;
+               d_paths[idpath*(r-1)+p-1].y = jc - jprec ;
+         }
+         jprec = jc ;
+       }
+       idpath++ ;
+  }
+  
+  // Q2
+  basepath += 6 ;
+  for (int a=0 ; a< 6 ; a++){         // les 6 angles 90,105,120,135,150,165
+       for (int p=0 ; p<r-1 ; p++){      // les r points
+         d_paths[idpath*(r-1)+p].x = -d_paths[(idpath - basepath)*(r-1)+p].y ;
+         d_paths[idpath*(r-1)+p].y =  d_paths[(idpath - basepath)*(r-1)+p].x ;
+         }
+       idpath++ ;
+  }
+
+  // Q3
+  basepath += 6 ;
+  for (int a=0 ; a< 6 ; a++){         // les 6 angles 180,195,210,225,240,255
+       for (int p=0 ; p<r-1 ; p++){      // les r points
+         d_paths[idpath*(r-1)+p].x = -d_paths[(idpath - basepath)*(r-1)+p].x ;
+         d_paths[idpath*(r-1)+p].y = -d_paths[(idpath - basepath)*(r-1)+p].y ;
+         }
+       idpath++ ;
+  }
+
+  // Q4
+  basepath += 6 ;
+  for (int a=0 ; a< 6 ; a++){         // les 6 angles 270,285,300,315,330,345
+       for (int p=0 ; p<r-1 ; p++){      // les r points
+         d_paths[idpath*(r-1)+p].x =  d_paths[(idpath - basepath)*(r-1)+p].y ;
+         d_paths[idpath*(r-1)+p].y = -d_paths[(idpath - basepath)*(r-1)+p].x ;
+         }
+       idpath++ ;
+  }
+}
 
 
+/**
+ *
+ * \brief calcule l'estimation initiale
+ * \author zulu - AND
+ *
+ * \param[in] L         Largeur de l'image
+ * \param[in] H         Hauteur de l'image
+ * \param[in] r         coté de la fenêtre de moyenneage
+ *
+ * \param[out] d_estim  Image estimee 0
+ *
+ * Version texture : l'img originale est supposée en texture.
+ * L'estimation réalisée correspond a un moyenneur de 'rayon' r
+ * Execution sur des blocs de threads 2D et une grille 2D
+ * selon les dimensions de l'image.
+ * 
+ */
+__global__ void kernel_neutre_img2estim(unsigned int *d_estim, unsigned int L, unsigned int H){
 
 
+  // coordonnes du point dans l'image
+  unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
+  unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
+  unsigned int pos = i*L +j ;
+
+  d_estim[ pos ] = tex2D(tex_img_in, j, i) ;
+  
+}
+
+
+/**
+ *
+ * \brief calcule l'estimation initiale
+ * \author zulu - AND
+ *
+ * \param[in] L         Largeur de l'image
+ * \param[in] H         Hauteur de l'image
+ * \param[in] r         coté de la fenêtre de moyenneage
+ *
+ * \param[out] d_estim  Image estimee 0
+ *
+ * Version texture : l'img originale est supposée en texture.
+ * L'estimation réalisée correspond a un moyenneur de 'rayon' r
+ * Execution sur des blocs de threads 2D et une grille 2D
+ * selon les dimensions de l'image.
+ * 
+ */
 __global__ void kernel_init_estim_from_img_in(unsigned int * d_estim, unsigned int L, unsigned int H, unsigned int r){
   // coordonnes du point dans l'image
   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 __global__ void kernel_init_estim_from_img_in(unsigned int * d_estim, unsigned int L, unsigned int H, unsigned int r){
   // coordonnes du point dans l'image
   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
@@ -91,11 +221,33 @@ __global__ void kernel_init_estim_from_img_in(unsigned int * d_estim, unsigned i
        d_estim[ pos ] = ng/((2*r+1)*(2*r+1)) ;
   }
   /*
        d_estim[ pos ] = ng/((2*r+1)*(2*r+1)) ;
   }
   /*
+       // pour les bords  : pas de traitement 
   else
   d_estim[ pos ] = tex2D(tex_img_in, j, i) ;
   */
 }
 
   else
   d_estim[ pos ] = tex2D(tex_img_in, j, i) ;
   */
 }
 
+
+
+/**
+ *
+ * \brief calcule l'estimation initiale
+ * \author zulu - AND
+ *
+ * \param[in] d_data    image originale
+ * \param[in] L         Largeur de l'image
+ * \param[in] H         Hauteur de l'image
+ * \param[in] r         coté de la fenêtre de moyenneage
+ *
+ * \param[out] d_estim  Image estimee 0
+ *
+ * Version global mem : l'img originale est en mémoire globale, passée en param.
+ * L'estimation réalisée correspond a un moyenneur de 'rayon' r
+ * Execution sur des blocs de threads 2D et une grille 2D
+ * selon les dimensions de l'image.
+ * Moins rapide que les 2 autres solutions.
+ * 
+ */
 __global__ void kernel_init_estim_from_img_in_global_mem(unsigned int * d_data, unsigned int * d_estim,
                                                                                                                 unsigned int L, unsigned int H, unsigned int r){
   // coordonnes du point dans l'image
 __global__ void kernel_init_estim_from_img_in_global_mem(unsigned int * d_data, unsigned int * d_estim,
                                                                                                                 unsigned int L, unsigned int H, unsigned int r){
   // coordonnes du point dans l'image
@@ -115,7 +267,28 @@ __global__ void kernel_init_estim_from_img_in_global_mem(unsigned int * d_data,
   } 
 }
 
   } 
 }
 
-__global__ void kernel_estim_next_step_global_mem(unsigned int * d_estim, unsigned int * d_lniv, unsigned int L, unsigned int H, unsigned int p){
+
+
+/**
+ *
+ * \brief calcule les niveaux de gris de l'estimation n+1
+ * \author zulu - AND
+ * 
+ * \param[in] L         Largeur de l'image
+ * \param[in] H         Hauteur de l'image
+ * \param[in] p         poids du terme lniv
+ *
+ * \param[out] d_estim  Image estimee n+1 
+ *
+ * Version mixte : l'img originale est supposee en texture,
+ * l'img lniv en mémoire globale, passée en param.
+ * Cela évite la copie en texture de l'img lniv à chaque itération.
+ * Execution sur des blocs de threads 2D et une grille 2D
+ * selon les dimensions de l'image.
+ * Moins rapide que 'texture' mais plus rapide que 'globalmem'
+ * 
+ */
+__global__ void kernel_estim_next_step_hybrid(unsigned int * d_estim, unsigned int * d_lniv, unsigned int L, unsigned int H, unsigned int p){
   // coordonnes du point dans l'image
   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
   // coordonnes du point dans l'image
   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
@@ -125,6 +298,23 @@ __global__ void kernel_estim_next_step_global_mem(unsigned int * d_estim, unsign
   
 }
 
   
 }
 
+/**
+ *
+ * \brief calcule les niveaux de gris de l'estimation n+1
+ * \author zulu - AND
+ * 
+ * \param[in] L         Largeur de l'image
+ * \param[in] H         Hauteur de l'image
+ * \param[in] p         poids du terme lniv
+ *
+ * \param[out] d_estim  Image estimee n+1 
+ *
+ * Version texture : Les donnees (img originale, img lniv) sont supposees en textures.
+ * Execution sur des blocs de threads 2D et une grille 2D
+ * selon les dimensions de l'image.
+ * Plus rapide que les 2 autres solutions
+ *
+ */
 __global__ void kernel_estim_next_step_texture(unsigned int * d_estim, unsigned int L, unsigned int H, unsigned int p){
   // coordonnes du point dans l'image
   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
 __global__ void kernel_estim_next_step_texture(unsigned int * d_estim, unsigned int L, unsigned int H, unsigned int p){
   // coordonnes du point dans l'image
   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
@@ -135,6 +325,24 @@ __global__ void kernel_estim_next_step_texture(unsigned int * d_estim, unsigned
   
 }
 
   
 }
 
+/**
+ *
+ * \brief calcule les niveaux de gris de l'estimation n+1
+ * \author zulu - AND
+ * 
+ * \param[in] d_lniv    Image des lniv n
+ * \param[in] d_data    Image originale
+ * \param[in] L         Largeur de l'image
+ * \param[in] H         Hauteur de l'image
+ * \param[in] p         poids du terme lniv
+ *
+ * \param[out] d_estim  Image estimee n+1 
+ *
+ * Version mémoire globale : les données sont passées en params.
+ * Execution sur des blocs de threads 2D et une grille 2D
+ * selon les dimensions de l'image.
+ * 
+ */
 __global__ void kernel_estim_next_step_global_mem(unsigned int * d_estim, unsigned int * d_lniv, unsigned int * d_data,
                                                                                                  unsigned int L, unsigned int H, unsigned int p){
   // coordonnes du point dans l'image
 __global__ void kernel_estim_next_step_global_mem(unsigned int * d_estim, unsigned int * d_lniv, unsigned int * d_data,
                                                                                                  unsigned int L, unsigned int H, unsigned int p){
   // coordonnes du point dans l'image
@@ -146,94 +354,57 @@ __global__ void kernel_estim_next_step_global_mem(unsigned int * d_estim, unsign
   
 }
 
   
 }
 
-__global__ void kernel_levelines_global_mem(unsigned int * img_in, unsigned int * img_out, unsigned int L, unsigned int H)
-{
-  // coordonnes du point dans l'image
-  unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
-  unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
-  //unsigned int spos = threadIdx.x * blockDim.y + threadIdx.y ;
-  
-  // nb de points par chemin
-  int lpath =  PSIZE_J ;
-  unsigned int ic, jc, zc, pos = i*L+j;
-  int idpath, idpix ;
-  unsigned int mse_min, mse_cur, val ;
-  uint2 mse ;
-  
-  
-  if((i>lpath)&&(i<H-lpath)&&(j>lpath)&&(j<L-lpath)){
-       for( idpath=0; idpath < PSIZE_I ; idpath++) {
-         ic = i ;
-         jc = j ;
-         pos = ic*L + jc ;
-         zc = img_in[ pos ] ;
-         mse.x = zc ;
-         mse.y = zc*zc ;
-         for( idpix=0; idpix < lpath-1 ; idpix++ ) {
-               ic += pathDi[idpath][idpix] ;
-               jc += pathDj[idpath][idpix] ;
-               pos = ic*L + jc ;
-               zc = img_in[ pos ] ;
-               mse.x += zc ;
-               mse.y += zc*zc ; 
-         }
-         // critere de selection du chemin ( SUM_(X2) - SUM_(X)2 / lpath )
-         // a ameliorer pour vitesse
-         mse_cur = ( mse.y - ( mse.x / lpath ) * mse.x ) ;
-         if (idpath == 0) {
-               mse_min = mse_cur ;
-               val = mse.x ;
-         } else {
-               if ( mse_cur < mse_min )  {
-                 mse_min = mse_cur ;
-                 val = mse.x ; 
-               }
-         } 
-       }
-       img_out[ i*L + j ] = val / lpath ; 
-  }
-}
-
-__global__ void kernel_levelines_texture(unsigned int * img_out, unsigned int L, unsigned int H)
+/**
+ *
+ * \brief determine les lniv en chaque point de l'image
+ * \author zulu - AND
+ *
+ * \param[in] L         Largeur de l'image
+ * \param[in] H         Hauteur de l'image
+ * \param[in] r         longueur des segments
+ *
+ * \param[out] img_out  image des lniv 
+ *
+ * Execution sur des blocs de threads 2D et une grille 2D
+ * selon les dimensions de l'image.
+ * L'image d'entrée doit être au préalable en mémoire texture pointée par "tex_img_estim".
+ * Les matrices des chemins sont, elles, pointées par "tex_paths"
+ * Cette version ne fournit pas les indices des chemins pour les tracé éventuel des lniv.
+ */
+__global__ void kernel_levelines_texture(unsigned int * img_out, unsigned int L, unsigned int H, unsigned int r)
 {
   // coordonnes du point dans l'image
   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
 {
   // coordonnes du point dans l'image
   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
-  //unsigned int spos = threadIdx.x * blockDim.y + threadIdx.y ;
-  
+    
   // nb de points par chemin
   // nb de points par chemin
-  int lpath =  PSIZE_J ;
-  unsigned int ic, jc, zc ;
+  int lpath =  r ;
+  unsigned int ic, jc, zc, z ;
   int idpath, idpix ;
   unsigned int mse_min, mse_cur, val ;
   uint2 mse ;
   
   
   int idpath, idpix ;
   unsigned int mse_min, mse_cur, val ;
   uint2 mse ;
   
   
-  if((i>lpath)&&(i<H-lpath)&&(j>lpath)&&(j<L-lpath)){
+  if((i>=lpath)&&(i<=H-lpath)&&(j>=lpath)&&(j<=L-lpath)){
+       z = tex2D(tex_img_estim, j, i) ;
        for( idpath=0; idpath < PSIZE_I ; idpath++) {
          ic = i ;
          jc = j ;
        for( idpath=0; idpath < PSIZE_I ; idpath++) {
          ic = i ;
          jc = j ;
-         zc = tex2D(tex_img_estim, j, i) ;
-         mse.x = zc ;
-         mse.y = zc*zc ;
+         mse.x = z ;
+         mse.y = z*z ;
          for( idpix=0; idpix < lpath-1 ; idpix++ ) {
          for( idpix=0; idpix < lpath-1 ; idpix++ ) {
-               ic += pathDi[idpath][idpix] ;
-               jc += pathDj[idpath][idpix] ;
+               ic += tex2D(tex_paths, idpix, idpath).x ;
+               jc += tex2D(tex_paths, idpix, idpath).y ;
                zc = tex2D(tex_img_estim, jc, ic) ;
                mse.x += zc ;
                mse.y += zc*zc ; 
          }
          // critere de selection du chemin ( SUM_(X2) - SUM_(X)2 / lpath )
                zc = tex2D(tex_img_estim, jc, ic) ;
                mse.x += zc ;
                mse.y += zc*zc ; 
          }
          // critere de selection du chemin ( SUM_(X2) - SUM_(X)2 / lpath )
-         // a ameliorer pour vitesse
+         // TODO cherchera  ameliorer pour vitesse
          mse_cur = ( mse.y - ( mse.x / lpath ) * mse.x ) ;
          mse_cur = ( mse.y - ( mse.x / lpath ) * mse.x ) ;
-         if (idpath == 0) {
+         if ( (idpath == 0) || (mse_cur < mse_min) ) {
                mse_min = mse_cur ;
                mse_min = mse_cur ;
-               val = mse.x ;
-         } else {
-               if ( mse_cur < mse_min )  {
-                 mse_min = mse_cur ;
-                 val = mse.x ; 
-               }
+               val = mse.x ; 
          } 
        }
        img_out[ i*L + j ] = val / lpath ; 
          } 
        }
        img_out[ i*L + j ] = val / lpath ; 
@@ -241,60 +412,101 @@ __global__ void kernel_levelines_texture(unsigned int * img_out, unsigned int L,
 }
 
 
 }
 
 
+/**
+ *
+ * \brief determine les lniv en chaque point de l'image
+ * \author zulu - AND
+ *
+ * \param[in] L         Largeur de l'image
+ * \param[in] H         Hauteur de l'image
+ * \param[in] r         longueur des segments
+ *
+ * \param[out] img_out  image des lniv 
+ *
+ * Execution sur des blocs de threads 2D et une grille 2D
+ * selon les dimensions de l'image.
+ * L'image d'entrée doit être au préalable en mémoire texture pointée par "tex_img_estim".
+ * Les matrices des chemins sont, elles, pointées par "tex_paths"
+ * Cette version ne fournit pas les indices des chemins pour les tracé éventuel des lniv.
+ * Cette version tente d'utiliser la shared memory pour compenser la baisse de perf due aux chemins
+ * paramétrables non constants.
+ */
 
 
-
-__global__ void kernel_levelines_only_texture(unsigned int * img_out, unsigned int L, unsigned int H)
+__global__ void kernel_levelines_texture_smem(unsigned int * img_out, unsigned int L, unsigned int H, unsigned int r)
 {
 {
+  // coordonnées du point dans le bloc
+  unsigned int iib = threadIdx.x ;
+  unsigned int jib = threadIdx.y ;
   // coordonnes du point dans l'image
   // coordonnes du point dans l'image
-  unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
-  unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;;
-  //unsigned int spos = threadIdx.x * blockDim.y + threadIdx.y ;
-  
+  unsigned int i = blockIdx.x*blockDim.x + iib ;
+  unsigned int j = blockIdx.y*blockDim.y + jib ;
+    
   // nb de points par chemin
   // nb de points par chemin
-  int lpath =  PSIZE_J ;
-  unsigned int ic, jc ;
+  int lpath =  r ;
+  int ic, jc ;
   int idpath, idpix ;
   int idpath, idpix ;
-  unsigned int mse_min, mse_cur ;
-  
-  //extern __shared__ uint2 mse[] ;
-  uint2 mse ;
+  unsigned int val, mse_cur, mse_min, z, zc ;
+  uint2 mse_data ;
+
+  //__shared__ unsigned int val_img[16*16] ;
+
+  //val_img[jib*16+iib] = tex2D(tex_img_estim, j, i) ;
   
   
-  if((i>lpath)&&(i<H-lpath)&&(j>lpath)&&(j<L-lpath)){
+  if((i>=lpath)&&(i<=H-lpath)&&(j>=lpath)&&(j<=L-lpath)){
+       z = tex2D(tex_img_estim, j, i) ;
        for( idpath=0; idpath < PSIZE_I ; idpath++) {
          ic = i ;
          jc = j ;
        for( idpath=0; idpath < PSIZE_I ; idpath++) {
          ic = i ;
          jc = j ;
-         mse.x = tex2D(tex_img_in, i, j) ;
-         mse.y = tex2D(tex_img_in, i, j)*tex2D(tex_img_in, i, j) ;
+         mse_data.x = z ;
+         mse_data.y = z*z ;
+         mse_min = mse_data.y - mse_data.x/lpath*mse_data.y ;
          for( idpix=0; idpix < lpath-1 ; idpix++ ) {
          for( idpix=0; idpix < lpath-1 ; idpix++ ) {
-               ic += pathDi[idpath][idpix] ;
-               jc += pathDj[idpath][idpix] ;
-               mse.x += tex2D( tex_img_in, ic, jc ) ;
-               mse.y += tex2D( tex_img_in, ic, jc ) * tex2D( tex_img_in, ic, jc ) ; 
+               ic += tex2D(tex_paths, idpix, idpath).x ;
+               jc += tex2D(tex_paths, idpix, idpath).y ;
+               zc = tex2D(tex_img_estim, jc, ic) ;
+               mse_data.x += zc ;
+               mse_data.y += zc*zc ; 
          }
          // critere de selection du chemin ( SUM_(X2) - SUM_(X)2 / lpath )
          }
          // critere de selection du chemin ( SUM_(X2) - SUM_(X)2 / lpath )
-         // a ameliorer pour vitesse
-         mse_cur = ( mse.y - ( mse.x / lpath ) * mse.x ) ;
-         if (idpath > 0) {
-               if ( mse_cur < mse_min )  {
-                 mse_min = mse_cur ;
-               }
-         } else {
+         // TODO cherchera  ameliorer pour vitesse
+         mse_cur = ( mse_data.y - ( mse_data.x / lpath ) * mse_data.x ) ;
+         if ( mse_cur < mse_min ){
                mse_min = mse_cur ;
                mse_min = mse_cur ;
-         }
+               val = mse_data.x ;
+         } 
        }
        }
-       img_out[ i*L + j ] = mse_min / lpath ; 
+       img_out[ i*L + j ] = val / lpath ; 
   }
 }
 
   }
 }
 
-
+/**
+ *
+ * \brief trace les segments sur un maillage carré
+ * \author zulu - AND
+ *
+ * \param[in] img_in    image d'entree
+ * \param[in] dir       tableaux des directions
+ * \param[in] L         Largeur de l'image
+ * \param[in] H         Hauteur de l'image
+ * \param[in] pas       coté du maillage
+ * \param[in] ng        niveau de gris des segments
+ * \param[in] r         longueur des segments
+ *
+ * \param[out] img_out  image + les segments superposés
+ *
+ * Kernel trivial. Ne trace rien sur les bords.
+ * execution sur des blocs de threads 2D et une grille 2D
+ * selon les dimensions de l'image
+ */
 __global__ void kernel_trace_levelines(unsigned int * img_in, unsigned int * dir, unsigned int * img_out,
 __global__ void kernel_trace_levelines(unsigned int * img_in, unsigned int * dir, unsigned int * img_out,
-                                                                          unsigned int L, unsigned int H, unsigned int pas, unsigned int ng){
+                                                                          unsigned int L, unsigned int H, unsigned int pas, unsigned int ng,
+                                                                          unsigned int r ){
   // coordonnes du point dans l'image
   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;;
     
   // nb de points par chemin
   // coordonnes du point dans l'image
   unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
   unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;;
     
   // nb de points par chemin
-  int lpath =  PSIZE_J ;
+  int lpath =  r ;
   unsigned int ic, jc, idpix ;
   unsigned int idpath  ;
 
   unsigned int ic, jc, idpix ;
   unsigned int idpath  ;
 
@@ -306,8 +518,8 @@ __global__ void kernel_trace_levelines(unsigned int * img_in, unsigned int * dir
        idpath = dir[ic*L+jc] ;
        img_out[ ic*L+jc ] = ng ;
        for ( idpix=0 ; idpix < lpath-1 ; idpix++ ){
        idpath = dir[ic*L+jc] ;
        img_out[ ic*L+jc ] = ng ;
        for ( idpix=0 ; idpix < lpath-1 ; idpix++ ){
-         ic += pathDi[idpath][idpix] ;
-         jc += pathDj[idpath][idpix] ;
+         ic += tex2D(tex_paths, idpix, idpath).x ; // pathDi[idpath][idpix] ;
+         jc += tex2D(tex_paths, idpix, idpath).y ; // pathDj[idpath][idpix] ;
          img_out[ ic*L + jc ] = ng ;
          }
   }
          img_out[ ic*L + jc ] = ng ;
          }
   }