#define PI 3.1415926

__constant__ float masque[256] ;
/* noyau gaussien 3x3 *//*
__constant__ float noyau[] = {1.0/16,2.0/16,1.0/16,
							2.0/16,4.0/16,2.0/16,
							1.0/16,2.0/16,1.0/16} ;
							__constant__ int rnoyau  = 1 ;*/
						
/* noyau gaussien 5x5 *//*
__constant__ int noyau[] = {1,4,6,4,1, 4,16,24,16,4, 6,24,36,24,6, 4,16,24,16,4, 1,4,6,4,1} ;
__constant__ int sumNoyau = 256 ;
__constant__ int rnoyau  = 2 ;
__constant__ int NNoyau = 25 ;
						*/
/* noyau 5x5*/
__constant__ float noyau[] = {1.0/25,1.0/25,1.0/25,1.0/25,1.0/25,
							  1.0/25,1.0/25,1.0/25,1.0/25,1.0/25,
							  1.0/25,1.0/25,1.0/25,1.0/25,1.0/25,
							  1.0/25,1.0/25,1.0/25,1.0/25,1.0/25,
							  1.0/25,1.0/25,1.0/25,1.0/25,1.0/25} ;
__constant__ int rnoyau  = 2 ;
			  
/* noyau 7x7 *//*
__device__ __constant__ float noyau[] = {1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,
										 1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,
										 1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,
										 1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,
										 1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,
										 1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,
										 1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49,1.0/49} ;
__device__ __constant__ int rnoyau  = 3 ;
			   */
//__device__ __constant__ int NNoyau = 49 ;

__device__ __constant__ int noyauH[] = {1,1,1,1,1,1,1} ;
__device__ __constant__ int sumNoyauH = 7 ;
__device__ __constant__ int rnoyauH  = 3 ;
// average 7x7
/*
__device__ __constant__ int noyauV[] = {1,1,1,1,1,1,1} ;
__device__ __constant__ int sumNoyauV = 7 ;
__device__ __constant__ int rnoyauV  = 3 ;
*/
// gaussian 5x5
__device__ __constant__ int noyauV[] = {1,4,6,4,1} ;
__device__ __constant__ int rnoyauV  = 2 ;
__constant__ int sumKV = 16 ;
//__device__ __constant__ int noyauV[] = {1,1,1} ;
//__constant__ int sumNoyauV = 16 ;

// declarations des textures 
texture<int, 2, cudaReadModeElementType> tex_img_in ;
texture<unsigned short, 2, cudaReadModeElementType> tex_img_ins ;
texture<unsigned char, 2, cudaReadModeElementType> tex_img_inc ;
texture<int, 2, cudaReadModeElementType> tex_img_estim ;
texture<float, 1, cudaReadModeElementType> tex_noyau ;
texture<float, 1, cudaReadModeElementType> tex_noyauConvof ;
texture<short, 1, cudaReadModeElementType> tex_noyauConvos ;
texture<int, 1, cudaReadModeElementType> tex_kern ;


/***************************************************************
 * kernel identite pour mesures de debit max effectif possible
 ***************************************************************/
__global__ void kernel_ident( unsigned short*output, int i_dim, int j_dim)
{  
  // coordonnees absolues du point de base
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  int idx = __mul24(i, j_dim) + j ;              // indice  dans l'image

  output[ idx ] = tex2D(tex_img_ins, j, i) ;

}

__global__ void kernel_ident2p( unsigned short  *output, int i_dim, int j_dim)
{  
  // coordonnees absolues du point de base
  int j = __umul24(__umul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ; 
  int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
   
  output[ __umul24(i, j_dim) + j    ] = tex2D(tex_img_inc, j, i) ;
  output[ __umul24(i, j_dim) + j +1 ] = tex2D(tex_img_inc, j+1, i) ;

}


__global__ void kernel_ident( unsigned char  *output, int i_dim, int j_dim)
{  
  // coordonnees absolues du point de base
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;

  output[ __umul24(i, j_dim) + j ] = tex2D(tex_img_inc, j, i) ;

}

__global__ void kernel_ident2p( unsigned char  *output, int i_dim, int j_dim)
{  
  // coordonnees absolues du point de base
  int j = __umul24(__umul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ; 
  int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
   
  output[ __umul24(i, j_dim) + j    ] = tex2D(tex_img_inc, j, i) ;
  output[ __umul24(i, j_dim) + j +1 ] = tex2D(tex_img_inc, j+1, i) ;

}

/**
 *
 * \brief calcule l'estimation initiale
 * \author zulu - AND
 *
 * \param[in] L         Largeur de l'image
 * \param[in] i_dim     Hauteur de l'image
 * \param[in] j_dim     coté de la fenêtre de moyenneage
 *
 * \param[out] output   Image estimee initiale
 *
 * 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_moyenne_shtex( unsigned int *output, unsigned int i_dim, unsigned int j_dim, unsigned int r)
{
  int idl, idc ;
  //coordonnées ds le bloc
  int ib = threadIdx.y ;
  int jb = threadIdx.x ;
  int idb_row = __mul24(ib,blockDim.x) + jb ;   // indice dans le bloc
  
  // coordonnees absolues du point
  int j = __mul24(blockIdx.x,blockDim.x) + jb ; 
  int i = __mul24(blockIdx.y,blockDim.y) + ib ;
  int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
  int off = __mul24(r,blockDim.x) ;             // offset du bloc dans le buffer en sh mem
  
  
  extern __shared__ unsigned int shmem[] ;
  volatile unsigned int * buffer = shmem ;
  volatile unsigned int * bufferV = &buffer[ blockDim.x*blockDim.y + 2*off ] ;
  
  
  /***********************************************************************************
   *               PASSE HORIZONTALE
   ***********************************************************************************/
  // charge les pixels de la bande au dessus du bloc
  if(ib<r)                                                           
	{
	  buffer[ idb_row ] = 0 ;
	  for (idc=j-r ; idc<=j+r ; idc++ )
		buffer[ idb_row  ] += tex2D(tex_img_in, idc, i-r) ;
	}

  // charge les pixels de la bande en dessous du bloc
  if (ib>=(blockDim.y-r))
	{
	  buffer[ idb_row + 2*off ] = 0 ;
	  for (idc=j-r ; idc<=j+r ; idc++ )
		buffer[ idb_row + 2*off ] += tex2D(tex_img_in, idc, i+r) ;
	}
  
  // charge les pixels du bloc
  buffer[ idb_row + off ] = 0 ;
  for (idc=j-r ; idc<=j+r ; idc++ )
	buffer[ idb_row + off ] += tex2D(tex_img_in, idc, i) ;
  
   __syncthreads() ;
  /*****************************************************************************
   *                             PASSE VERTICALE
   *****************************************************************************/
   bufferV[ idb_row ] = 0 ;
   for (idl=0 ; idl<=2*r ; idl++ )
	 bufferV[ idb_row ] += buffer[ (ib+idl)*blockDim.x + jb ] ;

  __syncthreads() ;
  
  output[idx] = bufferV[ idb_row ]/((2*r+1)*(2*r+1)) ; //pas juste sur les bords, mais c'est pô grave !
  
}

/**
 *
 * \brief calcule la convolution avec un noyau carre
 * \author zulu - AND
 *
 * \param[in] L         Largeur de l'image
 * \param[in] i_dim     Hauteur de l'image
 * \param[in] j_dim     coté de la fenêtre de moyenneage
 *
 * \param[out] output   Image estimee initiale
 *
 * Version texture : l'img originale est supposée en texture.
 * Execution sur des blocs de threads 2D et une grille 2D
 * correspondant aux dimensions de l'image.
 * 
 */

// convolution non separable
// image en texture et noyau en mem constante
__global__ void kernel_convoGene8( unsigned char  *output, int i_dim, int j_dim)
{
  int ic, jc ;
  int r = 2 ;
  int L=5 ;
  float outval0=0.0 ;
  
  // coordonnees absolues du point de base
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ; 
 
  for (ic=0 ; ic<L ; ic++)
	for(jc=0 ; jc<L ; jc++)
	  outval0 += masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic) ;

  
  output[ __mul24(i, j_dim) + j ] = outval0 ;
}

// convolution non separable 2 pixels/thread
// image en texture et noyau en mem constante
__global__ void kernel_convoGene8x2p( unsigned char  *output, int i_dim, int j_dim )
{
  int ic, jc ;
  int r = 2 ;
  int L=5 ;
  float outval0=0.0, outval1=0.0, outval2=0.0 ;
  
  // coordonnees absolues du point de base
  int j = __mul24( __mul24(blockIdx.x,blockDim.x) + threadIdx.x, 2) ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ; 
 
  for (ic=0 ; ic<L ; ic++)
	{
	  for(jc=1 ; jc<L-1 ; jc++)
		{
		  outval0 += masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic) ;		
		}
	}
  outval1 = outval0 ;
  outval2 = outval0 ;
  for (ic=0 ; ic<L ; ic++)
	{
	  outval1 += masque[ __mul24(ic,L) ]*tex2D(tex_img_inc, j-r, i-r+ic) ;
	  outval2 += masque[ __mul24(ic,L)+L-1 ]*tex2D(tex_img_inc, j-r+L, i-r+ic) ;		
	}
  
  output[ __mul24(i, j_dim) + j ] = outval1 ;
  output[ __mul24(i, j_dim) + j+1 ] = outval2 ;
}

// convolution non separable 2 pixels/thread
// image en texture et noyau en mem constante
__global__ void kernel_convoGene8rx2p( unsigned char  *output, int j_dim, int r )
{
  int ic, jc ;
  int L=2*r+1 ;
  //unsigned char pix ;
  float outval0=0.0, outval1=0.0 ;
  
  // coordonnees absolues du point de base
  int j = __umul24( __umul24(blockIdx.x,blockDim.x) + threadIdx.x, 2) ; 
  int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ; 
  int idx=__umul24(i, j_dim) + j ;
  
  for (ic=0 ; ic<L ; ic++)
	{
	  outval0 += masque[ __umul24(ic,L) ]*tex2D(tex_img_inc, j-r, i-r+ic) ; //bande gauche
	  //zone commune
	  for(jc=1 ; jc<L ; jc++)
		{
		  //pix = tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
		  outval0 += masque[ __umul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
		  outval1 += masque[ __umul24(ic,L)+jc-1 ]*tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
		}
	  outval1 += masque[ __umul24(ic,L)+L-1 ]*tex2D(tex_img_inc, j-r+L, i-r+ic) ; // bande droite      
	}
  output[ idx   ] = outval0 ;
  output[ idx+1 ] = outval1 ;
}

// convolution non separable 4 pixels/thread
// image en texture et noyau en mem constante
__global__ void kernel_convoGene8rx4p( unsigned char  *output, int j_dim, int r )
{
  int ic, jc ;
  int L=2*r+1 ;
  unsigned char pix ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
  
  // coordonnees absolues du point de base en haut a gauche
  int j = __mul24( __umul24( blockIdx.x, blockDim.x) + threadIdx.x, 2) ; 
  int i = __mul24( __umul24( blockIdx.y, blockDim.y) + threadIdx.y, 2) ; 

  //zone commune
  //forme generique contrib pixel => masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic)
  for (ic=1 ; ic<L ; ic++)
	{
	  for(jc=1 ; jc<L ; jc++)
		{
		  pix = tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
		  outval0 += masque[ __umul24(ic,L)  +jc   ]*pix ;
		  outval1 += masque[ __umul24(ic,L)  +jc-1 ]*pix ;
		  outval2 += masque[ __umul24(ic-1,L)+jc   ]*pix ;
		  outval3 += masque[ __umul24(ic-1,L)+jc-1 ]*pix ;
		}
	}
  // les blocs peripheriques verticaux
  for (ic=1 ; ic<L ; ic++)
	{
	  outval0 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r, i-r+ic) ; // bande A pour b0
	  outval2 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r, i-r+ic) ; // bande A pour b2
	  
	  outval1 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+1, i-r+ic) ; //bande C pour b1
	  outval3 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+1, i-r+ic) ; //bande C pour b3
	}
  // les blocs peripheriques horizontaux
  for (jc=1 ; jc<L ; jc++)
	{
	  outval0 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc, i-r) ; // bande B pour b0
	  outval1 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc, i-r) ; // bande B pour b1
	  
	  outval2 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc, i+r+1) ; //bande D pour b2
	  outval3 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc, i+r+1) ; //bande D pour b3
	}

  //les coins
  {
	outval0 += masque[ 0                   ]*tex2D(tex_img_inc, j-r  , i-r) ;   // haut gauche pour  b0
	outval1 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+1, i-r) ;   // haut droit pour b1
	outval2 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r  , i+r+1) ; // bas gauche pour b2
	outval3 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+1, i+r+1) ; //bas droit pour b3
  }

  
  // les 2 pixels hauts 
  output[ __umul24(i, j_dim) + j   ] = outval0 ;
  output[ __umul24(i, j_dim) + j+1 ] = outval1 ;
  // les 2 pixels bas
  output[ __umul24(i+1, j_dim) + j   ] = outval2 ;
  output[ __umul24(i+1, j_dim) + j+1 ] = outval3 ; 
}

// convolution non separable 4 pixels/thread
// image en texture et noyau en mem constante
__global__ void kernel_convoGene8rx8p( unsigned char  *output, int j_dim, int r )
{
  int ic, jc ;
  int L=2*r+1 ;
  unsigned char pix ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
  float outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  
  // coordonnees absolues du point de base en haut a gauche
  int j = ( __umul24( blockIdx.x, blockDim.x) + threadIdx.x)<< 2 ; 
  int i = ( __umul24( blockIdx.y, blockDim.y) + threadIdx.y)<< 1 ; 

  //zone commune
  //forme generique contrib pixel => masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic)
  for (ic=1 ; ic<L ; ic++)
	{
	  for(jc=1 ; jc<L ; jc++)
		{
		  pix = tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
		  outval0 += masque[ __umul24(ic,L)  +jc   ]*pix ;
		  outval1 += masque[ __umul24(ic,L)  +jc-1 ]*pix ;
		  outval2 += masque[ __umul24(ic-1,L)+jc   ]*pix ;
		  outval3 += masque[ __umul24(ic-1,L)+jc-1 ]*pix ;

		  pix = tex2D(tex_img_inc, j-r+jc+2, i-r+ic) ;
		  outval4 += masque[ __umul24(ic,L)  +jc   ]*pix ;
		  outval5 += masque[ __umul24(ic,L)  +jc-1 ]*pix ;
		  outval6 += masque[ __umul24(ic-1,L)+jc   ]*pix ;
		  outval7 += masque[ __umul24(ic-1,L)+jc-1 ]*pix ;
		  
		}
	}
  // les blocs peripheriques verticaux
  for (ic=1 ; ic<L ; ic++)
	{
	  outval0 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r, i-r+ic) ; // bande A pour b0
	  outval2 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r, i-r+ic) ; // bande A pour b2
	  
	  outval1 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+1, i-r+ic) ; //bande C pour b1
	  outval3 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+1, i-r+ic) ; //bande C pour b3

	  outval4 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r+2, i-r+ic) ; // bande A pour b0
	  outval5 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r+2, i-r+ic) ; // bande A pour b2
	  
	  outval6 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+3, i-r+ic) ; //bande C pour b1
	  outval7 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+3, i-r+ic) ; //bande C pour b3
	}
  // les blocs peripheriques horizontaux
  for (jc=1 ; jc<L ; jc++)
	{
	  outval0 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc, i-r) ; // bande B pour b0
	  outval1 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc, i-r) ; // bande B pour b1
	  
	  outval2 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc, i+r+1) ; //bande D pour b2
	  outval3 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc, i+r+1) ; //bande D pour b3

	  outval4 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc+2, i-r) ; // bande B pour b0
	  outval5 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc+2, i-r) ; // bande B pour b1
	  
	  outval6 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc+2, i+r+1) ; //bande D pour b2
	  outval7 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc+2, i+r+1) ; //bande D pour b3
	}

  //les coins
  {
	outval0 += masque[ 0                   ]*tex2D(tex_img_inc, j-r  , i-r) ;   // haut gauche pour  b0
	outval1 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+1, i-r) ;   // haut droit pour b1
	outval2 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r  , i+r+1) ; // bas gauche pour b2
	outval3 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+1, i+r+1) ; //bas droit pour b3

	outval4 += masque[ 0                   ]*tex2D(tex_img_inc, j-r+2, i-r) ;   // haut gauche pour  b0
	outval5 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+3, i-r) ;   // haut droit pour b1
	outval6 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r+2, i+r+1) ; // bas gauche pour b2
	outval7 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+3, i+r+1) ; //bas droit pour b3
  }

  
  // les 2 pixels hauts 
  output[ __umul24(i, j_dim) + j   ] = outval0 ;
  output[ __umul24(i, j_dim) + j+1 ] = outval1 ;
  output[ __umul24(i, j_dim) + j+2 ] = outval4 ;
  output[ __umul24(i, j_dim) + j+3 ] = outval5 ;
  // les 2 pixels bas
  output[ __umul24(i+1, j_dim) + j   ] = outval2 ;
  output[ __umul24(i+1, j_dim) + j+1 ] = outval3 ;
  output[ __umul24(i+1, j_dim) + j+2 ] = outval6 ;
  output[ __umul24(i+1, j_dim) + j+3 ] = outval7 ;
   
}

// convolution non separable 8 pixels/thread
// image en texture et noyau en mem constante
__global__ void kernel_convoGene8x8pL3( unsigned char  *output, int j_dim )
{
  int ic, base ;
  int L=3 ;
  unsigned char pix ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
  float outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  
  // coordonnees absolues du point de base en haut a gauche
  int j = ( __umul24( blockIdx.x, blockDim.x) + threadIdx.x)<< 3 ; 
  int i = ( __umul24( blockIdx.y, blockDim.y) + threadIdx.y) ; 

  //zone commune
  //forme generique contrib pixel => masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic)
  for (ic=0 ; ic<L ; ic++)
	{
	  base = __umul24(ic,L) ;
	  pix = tex2D(tex_img_inc, j+1, i-1+ic) ;
	  outval0 += masque[ base +2 ]*pix ;
	  outval1 += masque[ base +1 ]*pix ;
	  outval2 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j+2, i-1+ic) ;
	  outval1 += masque[ base +2 ]*pix ;
	  outval2 += masque[ base +1 ]*pix ;
	  outval3 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j+3, i-1+ic) ;
	  outval2 += masque[ base +2 ]*pix ;
	  outval3 += masque[ base +1 ]*pix ;
	  outval4 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j+4, i-1+ic) ;
	  outval3 += masque[ base +2 ]*pix ;
	  outval4 += masque[ base +1 ]*pix ;
	  outval5 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j+5, i-1+ic) ;
	  outval4 += masque[ base +2 ]*pix ;
	  outval5 += masque[ base +1 ]*pix ;
	  outval6 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j+6, i-1+ic) ;
	  outval5 += masque[ base +2 ]*pix ;
	  outval6 += masque[ base +1 ]*pix ;
	  outval7 += masque[ base    ]*pix ;
  
	  pix = tex2D(tex_img_inc, j, i-1+ic) ;
	  outval0 += masque[ base +1 ]*pix ;
	  outval1 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j-1, i-1+ic) ;
	  outval0 += masque[ base  ]*pix ;

	  pix = tex2D(tex_img_inc, j+7, i-1+ic) ;
	  outval6 += masque[ base +2 ]*pix ;
	  outval7 += masque[ base +1 ]*pix ;
	  pix = tex2D(tex_img_inc, j+8, i-1+ic) ;
	  outval7 += masque[ base +2 ]*pix ;
	}
  
  // les 2 pixels hauts
  base = __umul24(i, j_dim) + j ;
  output[ base++ ] = outval0 ;
  output[ base++ ] = outval1 ;
  output[ base++ ] = outval2 ;
  output[ base++ ] = outval3;
  output[ base++ ] = outval4;
  output[ base++ ] = outval5;
  output[ base++ ] = outval6 ;
  output[ base++ ] = outval7 ;
   
}





// convolution non separable 8 pixels/thread
// image en texture et noyau en mem constante
__global__ void kernel_convoGene8x8pL5( unsigned char  *output, int j_dim )
{
  int ic, base ;
  int L=5 ;
  unsigned char pix ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
  float outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  
  // coordonnees absolues du point de base en haut a gauche
  int j = ( __umul24( blockIdx.x, blockDim.x) + threadIdx.x)<< 3 ; 
  int i = ( __umul24( blockIdx.y, blockDim.y) + threadIdx.y) ; 

  //zone commune
  //forme generique contrib pixel => masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic)
  for (ic=0 ; ic<L ; ic++)
	{
	  base = __umul24(ic,L) ;
	  pix = tex2D(tex_img_inc, j+2, i-2+ic) ;
	  outval0 += masque[ base +4 ]*pix ;
	  outval1 += masque[ base +3 ]*pix ;
	  outval2 += masque[ base +2 ]*pix ;
	  outval3 += masque[ base +1 ]*pix ;
	  outval4 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j+3, i-2+ic) ;
	  outval1 += masque[ base +4 ]*pix ;
	  outval2 += masque[ base +3 ]*pix ;
	  outval3 += masque[ base +2 ]*pix ;
	  outval4 += masque[ base +1 ]*pix ;
	  outval5 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j+4, i-2+ic) ;
	  outval2 += masque[ base +4 ]*pix ;
	  outval3 += masque[ base +3 ]*pix ;
	  outval4 += masque[ base +2 ]*pix ;
	  outval5 += masque[ base +1 ]*pix ;
	  outval6 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j+5, i-2+ic) ;
	  outval3 += masque[ base +4 ]*pix ;
	  outval4 += masque[ base +3 ]*pix ;
	  outval5 += masque[ base +2 ]*pix ;
	  outval6 += masque[ base +1 ]*pix ;
	  outval7 += masque[ base    ]*pix ;

	  pix = tex2D(tex_img_inc, j+1, i-2+ic) ;
	  outval0 += masque[ base +3 ]*pix ;
	  outval1 += masque[ base +2 ]*pix ;
	  outval2 += masque[ base +1 ]*pix ;
	  outval3 += masque[ base    ]*pix ;
	  
	  pix = tex2D(tex_img_inc, j, i-2+ic) ;
	  outval0 += masque[ base +2 ]*pix ;
	  outval1 += masque[ base +1 ]*pix ;
	  outval2 += masque[ base    ]*pix ;

	  pix = tex2D(tex_img_inc, j-1, i-2+ic) ;
	  outval0 += masque[ base +1 ]*pix ;
	  outval1 += masque[ base    ]*pix ;

	  pix = tex2D(tex_img_inc, j-2, i-2+ic) ;
	  outval0 += masque[ base    ]*pix ;

	  pix = tex2D(tex_img_inc, j+6, i-2+ic) ;
	  outval7 += masque[ base +1 ]*pix ;
	  outval6 += masque[ base +2 ]*pix ;
	  outval5 += masque[ base +3 ]*pix ;
	  outval4 += masque[ base +4 ]*pix ;
	  
	  pix = tex2D(tex_img_inc, j+7, i-2+ic) ;
	  outval7 += masque[ base +2 ]*pix ;
	  outval6 += masque[ base +3 ]*pix ;
	  outval5 += masque[ base +4 ]*pix ;

	  pix = tex2D(tex_img_inc, j+8, i-2+ic) ;
	  outval7 += masque[ base +3 ]*pix ;
	  outval6 += masque[ base +4 ]*pix ;

	  pix = tex2D(tex_img_inc, j+9, i-2+ic) ;
	  outval7 += masque[ base +4 ]*pix ;
	}
  
  // les 8 pixels
  base = __umul24(i, j_dim) + j ;
  output[ base++ ] = outval0 ;
  output[ base++ ] = outval1 ;
  output[ base++ ] = outval2 ;
  output[ base++ ] = outval3;
  output[ base++ ] = outval4;
  output[ base++ ] = outval5;
  output[ base++ ] = outval6 ;
  output[ base++ ] = outval7 ;
   
}

__global__ void kernel_convoGene8x8pL7( unsigned char  *output, int j_dim )
{
  int ic, base ;
  int L=7 ;
  unsigned char pix ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
  float outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  
  // coordonnees absolues du point de base en haut a gauche
  int j = ( __umul24( blockIdx.x, blockDim.x) + threadIdx.x)<< 3 ; 
  int i = ( __umul24( blockIdx.y, blockDim.y) + threadIdx.y) ; 

  //zone commune
  //forme generique contrib pixel => masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic)
  for (ic=0 ; ic<L ; ic++)
	{
	  base = __umul24(ic,L) ;
	  pix = tex2D(tex_img_inc, j+3, i-3+ic) ;
	  outval0 += masque[ base +6 ]*pix ;
	  outval1 += masque[ base +5 ]*pix ;
	  outval2 += masque[ base +4 ]*pix ;
	  outval3 += masque[ base +3 ]*pix ;
	  outval4 += masque[ base +2 ]*pix ;
	  outval5 += masque[ base +1 ]*pix ;
	  outval6 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j+4, i-3+ic) ;
	  outval1 += masque[ base +6 ]*pix ;
	  outval2 += masque[ base +5 ]*pix ;
	  outval3 += masque[ base +4 ]*pix ;
	  outval4 += masque[ base +3 ]*pix ;
	  outval5 += masque[ base +2 ]*pix ;
	  outval6 += masque[ base +1 ]*pix ;
	  outval7 += masque[ base    ]*pix ;
	  
	  pix = tex2D(tex_img_inc, j+2, i-3+ic) ;
	  outval0 += masque[ base +5 ]*pix ;
	  outval1 += masque[ base +4 ]*pix ;
	  outval2 += masque[ base +3 ]*pix ;
	  outval3 += masque[ base +2 ]*pix ;
	  outval4 += masque[ base +1 ]*pix ;
	  outval5 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j+1, i-3+ic) ;
	  outval0 += masque[ base +4 ]*pix ;
	  outval1 += masque[ base +3 ]*pix ;
	  outval2 += masque[ base +2 ]*pix ;
	  outval3 += masque[ base +1 ]*pix ;
	  outval4 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j, i-3+ic) ;
	  outval0 += masque[ base +3 ]*pix ;
	  outval1 += masque[ base +2 ]*pix ;
	  outval2 += masque[ base +1 ]*pix ;
	  outval3 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j-1, i-3+ic) ;
	  outval0 += masque[ base +2 ]*pix ;
	  outval1 += masque[ base +1 ]*pix ;
	  outval2 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j-2, i-3+ic) ;
	  outval0 += masque[ base +1 ]*pix ;
	  outval1 += masque[ base    ]*pix ;
	  pix = tex2D(tex_img_inc, j-3, i-3+ic) ;
	  outval0 += masque[ base ]*pix ;
	  

	  pix = tex2D(tex_img_inc, j+5, i-3+ic) ;
	  outval7 += masque[ base +1 ]*pix ;
	  outval6 += masque[ base +2 ]*pix ;
	  outval5 += masque[ base +3 ]*pix ;
	  outval4 += masque[ base +4 ]*pix ;
	  outval3 += masque[ base +5 ]*pix ;
	  outval2 += masque[ base +6 ]*pix ;
	  pix = tex2D(tex_img_inc, j+6, i-3+ic) ;
	  outval7 += masque[ base +2 ]*pix ;
	  outval6 += masque[ base +3 ]*pix ;
	  outval5 += masque[ base +4 ]*pix ;
	  outval4 += masque[ base +5 ]*pix ;
	  outval3 += masque[ base +6 ]*pix ;
	  pix = tex2D(tex_img_inc, j+7, i-3+ic) ;
	  outval7 += masque[ base +3 ]*pix ;
	  outval6 += masque[ base +4 ]*pix ;
	  outval5 += masque[ base +5 ]*pix ;
	  outval4 += masque[ base +6 ]*pix ;
	  pix = tex2D(tex_img_inc, j+8, i-3+ic) ;
	  outval7 += masque[ base +4 ]*pix ;
	  outval6 += masque[ base +5 ]*pix ;
	  outval5 += masque[ base +6 ]*pix ;
	  pix = tex2D(tex_img_inc, j+9, i-3+ic) ;
	  outval7 += masque[ base +5 ]*pix ;
	  outval6 += masque[ base +6 ]*pix ;
	  pix = tex2D(tex_img_inc, j+10, i-3+ic) ;
	  outval7 += masque[ base +6 ]*pix ;
	  
	}
  
  // les 8 pixels
  base = __umul24(i, j_dim) + j ;
  output[ base++ ] = outval0 ;
  output[ base++ ] = outval1 ;
  output[ base++ ] = outval2 ;
  output[ base++ ] = outval3;
  output[ base++ ] = outval4;
  output[ base++ ] = outval5;
  output[ base++ ] = outval6 ;
  output[ base++ ] = outval7 ; 
}



// convolution non separable 4 pixels/thread
// image en texture et noyau en mem constante
__global__ void kernel_convoGene8rx16p( unsigned char  *output, int j_dim, int r )
{
  int ic, jc ;
  int L=2*r+1, baseIdx ;
  unsigned char pix ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
  float outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  float outval8=0.0, outval9=0.0, outval10=0.0, outval11=0.0 ;
  float outval12=0.0, outval13=0.0, outval14=0.0, outval15=0.0 ;
  
  // coordonnees absolues du point de base en haut a gauche
  int j = ( __umul24( blockIdx.x, blockDim.x) + threadIdx.x)<< 3 ; 
  int i = ( __umul24( blockIdx.y, blockDim.y) + threadIdx.y)<< 1 ; 

  //zone commune
  //forme generique contrib pixel => masque[ __mul24(ic,L)+jc ]*tex2D(tex_img_inc, j-r+jc, i-r+ic)
  for (ic=1 ; ic<L ; ic++)
	{
	  for(jc=1 ; jc<L ; jc++)
		{
		  pix = tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
		  baseIdx = __umul24(ic-1,L)+jc-1 ;
		  outval0 += masque[ baseIdx+L+1 ]*pix ;
		  outval1 += masque[ baseIdx+L   ]*pix ;
		  outval2 += masque[ baseIdx+1   ]*pix ;
		  outval3 += masque[ baseIdx     ]*pix ;

		  pix = tex2D(tex_img_inc, j-r+jc+2, i-r+ic) ;
		  outval4 += masque[ baseIdx+L+1  ]*pix ;
		  outval5 += masque[ baseIdx+L    ]*pix ;
		  outval6 += masque[ baseIdx+1    ]*pix ;
		  outval7 += masque[ baseIdx      ]*pix ;

		  pix = tex2D(tex_img_inc, j-r+jc+4, i-r+ic) ;
		  outval8 += masque[ baseIdx+L+1 ]*pix ;
		  outval9 += masque[ baseIdx+L   ]*pix ;
		  outval10+= masque[ baseIdx+1   ]*pix ;
		  outval11+= masque[ baseIdx     ]*pix ;

		  pix = tex2D(tex_img_inc, j-r+jc+6, i-r+ic) ;
		  outval12 += masque[ baseIdx+L+1 ]*pix ;
		  outval13 += masque[ baseIdx+L   ]*pix ;
		  outval14 += masque[ baseIdx+1   ]*pix ;
		  outval15 += masque[ baseIdx     ]*pix ;
		  
		}
	}
  // les blocs peripheriques verticaux
  for (ic=1 ; ic<L ; ic++)
	{
	  outval0 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r, i-r+ic) ; // bande A pour b0
	  outval2 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r, i-r+ic) ; // bande A pour b2
	  
	  outval1 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+1, i-r+ic) ; //bande C pour b1
	  outval3 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+1, i-r+ic) ; //bande C pour b3

	  outval4 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r+2, i-r+ic) ; // bande A pour b0
	  outval6 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r+2, i-r+ic) ; // bande A pour b2
	  
	  outval5 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+3, i-r+ic) ; //bande C pour b1
	  outval7 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+3, i-r+ic) ; //bande C pour b3

	  outval8 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r+4, i-r+ic) ; // bande A pour b0
	  outval10 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r+4, i-r+ic) ; // bande A pour b2
	  
	  outval9 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+5, i-r+ic) ; //bande C pour b1
	  outval11 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+5, i-r+ic) ; //bande C pour b3

	  outval12 += masque[ __umul24(ic,L)   ]*tex2D(tex_img_inc, j-r+6, i-r+ic) ; // bande A pour b0
	  outval14 += masque[ __umul24(ic-1,L) ]*tex2D(tex_img_inc, j-r+6, i-r+ic) ; // bande A pour b2
	  
	  outval13 += masque[ __umul24(ic,L)  +L-1 ]*tex2D(tex_img_inc, j+r+7, i-r+ic) ; //bande C pour b1
	  outval15 += masque[ __umul24(ic-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+7, i-r+ic) ; //bande C pour b3
	}
  // les blocs peripheriques horizontaux
  for (jc=1 ; jc<L ; jc++)
	{
	  outval0 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc, i-r) ; // bande B pour b0
	  outval1 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc, i-r) ; // bande B pour b1
	  
	  outval2 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc, i+r+1) ; //bande D pour b2
	  outval3 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc, i+r+1) ; //bande D pour b3

	  outval4 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc+2, i-r) ; // bande B pour b0
	  outval5 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc+2, i-r) ; // bande B pour b1
	  
	  outval6 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc+2, i+r+1) ; //bande D pour b2
	  outval7 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc+2, i+r+1) ; //bande D pour b3

	  outval8 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc+4, i-r) ; // bande B pour b0
	  outval9 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc+4, i-r) ; // bande B pour b1
	  
	  outval10 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc+4, i+r+1) ; //bande D pour b2
	  outval11 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc+4, i+r+1) ; //bande D pour b3

	  outval12 += masque[ jc   ]*tex2D(tex_img_inc, j-r+jc+6, i-r) ; // bande B pour b0
	  outval13 += masque[ jc-1 ]*tex2D(tex_img_inc, j-r+jc+6, i-r) ; // bande B pour b1
	  
	  outval14 += masque[ __umul24(L-1,L) +jc   ]*tex2D(tex_img_inc, j-r+jc+6, i+r+1) ; //bande D pour b2
	  outval15 += masque[ __umul24(L-1,L) +jc-1 ]*tex2D(tex_img_inc, j-r+jc+6, i+r+1) ; //bande D pour b3
	}

  //les coins
  {
	outval0 += masque[ 0                   ]*tex2D(tex_img_inc, j-r  , i-r) ;   // haut gauche pour  b0
	outval1 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+1, i-r) ;   // haut droit pour b1
	outval2 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r  , i+r+1) ; // bas gauche pour b2
	outval3 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+1, i+r+1) ; //bas droit pour b3

	outval4 += masque[ 0                   ]*tex2D(tex_img_inc, j-r+2, i-r) ;   // haut gauche pour  b0
	outval5 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+3, i-r) ;   // haut droit pour b1
	outval6 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r+2, i+r+1) ; // bas gauche pour b2
	outval7 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+3, i+r+1) ; //bas droit pour b3

	outval8 += masque[ 0                   ]*tex2D(tex_img_inc, j-r+4, i-r) ;   // haut gauche pour  b0
	outval9 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+5, i-r) ;   // haut droit pour b1
	outval10 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r+4, i+r+1) ; // bas gauche pour b2
	outval11 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+5, i+r+1) ; //bas droit pour b3

	outval12 += masque[ 0                   ]*tex2D(tex_img_inc, j-r+6, i-r) ;   // haut gauche pour  b0
	outval13 += masque[ L-1                 ]*tex2D(tex_img_inc, j+r+7, i-r) ;   // haut droit pour b1
	outval14 += masque[ __umul24(L-1,L)     ]*tex2D(tex_img_inc, j-r+6, i+r+1) ; // bas gauche pour b2
	outval15 += masque[ __umul24(L-1,L)+L-1 ]*tex2D(tex_img_inc, j+r+7, i+r+1) ; //bas droit pour b3
  }

  
  // les 2 pixels hauts
  int base = __umul24(i, j_dim) +j ; 
  output[ base++   ] = outval0 ;
  output[ base++ ] = outval1 ;
  output[ base++ ] = outval4 ;
  output[ base++ ] = outval5 ;
  output[ base++ ] = outval8 ;
  output[ base++ ] = outval9 ;
  output[ base++ ] = outval12 ;
  output[ base++ ] = outval13 ;
  // les 2 pixels bas
  base = __umul24(i+1, j_dim) +j ;
  output[ base++ ] = outval2 ;
  output[ base++ ] = outval3 ;
  output[ base++ ] = outval6 ;
  output[ base++ ] = outval7 ;
  output[ base++ ] = outval10 ;
  output[ base++ ] = outval11 ;
  output[ base++ ] = outval14 ;
  output[ base++ ] = outval15 ;

}



// convolution non separable
// image en texture et noyau en mem constante
__global__ void kernel_convoGene8r( unsigned char  *output, int i_dim, int j_dim, int r)
{
  int ic, jc ;
  int L=2*r+1 ;
  float outval0=0.0 ;
  
  // coordonnees absolues du point de base
  int j = __mul24( blockIdx.x,blockDim.x ) + threadIdx.x ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ; 
 
  for (ic=-r ; ic<=r ; ic++)
	for(jc=-r ; jc<=r ; jc++)
	  outval0 += masque[ __umul24(ic,L)+jc+r ]*tex2D(tex_img_inc, j+jc, i+ic) ;

  
  output[ __umul24(i, j_dim) + j ] = outval0 ;
}

// convolution non separable
// image en texture et noyau en mem constante
__global__ void kernel_convoGeneTex8r( unsigned char  *output, int i_dim, int j_dim, int r)
{
  int ic, jc ;
  int L=2*r+1 ;
  float outval0=0.0 ;
  
  // coordonnees absolues du point de base
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ; 
 
  for (ic=0 ; ic<L ; ic++)
	for(jc=0 ; jc<L ; jc++)
	  outval0 += tex1D(tex_noyau, __mul24(ic,L)+jc )*tex2D(tex_img_inc, j-r+jc, i-r+ic) ;
  
  output[ __mul24(i, j_dim) + j ] = outval0 ;
}


__global__ void kernel_convoGene3Reg8( unsigned char  *output, int i_dim, int j_dim )
{
  float outval0=0.0 ;
  float n0,n1,n2,n3,n4,n5,n6,n7,n8 ; 
  
  n0 = (1.0/45) ;
  n1 = (2.0/45) ;
  n2 = (3.0/45) ;
  n3 = (4.0/45) ;
  n4 = (5.0/45) ;
  n5 = (6.0/45) ;
  n6 = (7.0/45) ;
  n7 = (8.0/45) ;
  n8 = (9.0/45) ;
  
  // coordonnees absolues du point de base
  int j = __mul24(blockIdx.x, blockDim.x) + threadIdx.x ; 
  int i = __mul24(blockIdx.y, blockDim.y) + threadIdx.y ; 
    
  outval0 = n8*tex2D(tex_img_inc, j-1, i-1 ) 
	      + n7*tex2D(tex_img_inc, j  , i-1 ) 
	      + n6*tex2D(tex_img_inc, j+1, i-1 ) 
	      + n5*tex2D(tex_img_inc, j-1, i   ) 
	      + n4*tex2D(tex_img_inc, j  , i   ) 
          + n3*tex2D(tex_img_inc, j+1, i   ) 
	      + n2*tex2D(tex_img_inc, j-1, i+1 ) 
	      + n1*tex2D(tex_img_inc, j  , i+1 ) 
	      + n0*tex2D(tex_img_inc, j+1, i+1 ) ;
	      
  output[ __mul24(i, j_dim) + j  ] = (unsigned char) outval0 ;
}


// convolution non separable
// image en texture et noyau en registres
__global__ void kernel_convoGene5Reg8( unsigned char  *output, int i_dim, int j_dim)
{
  float outval0=0.0 ;
  float n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n13,n14,n15,n16,n17,n18,n19,n20,n21,n22,n23,n24 ; 
  
  n0 = (1.0/25) ;
  n1 = (1.0/25) ;
  n2 = (1.0/25) ;
  n3 = (1.0/25) ;
  n4 = (1.0/25) ;
  n5 = (1.0/25) ;
  n6 = (1.0/25) ;
  n7 = (1.0/25) ;
  n8 = (1.0/25) ;
  n9 = (1.0/25) ;
  n10 = (1.0/25) ;
  n11 = (1.0/25) ;
  n12 = (1.0/25) ;
  n13 = (1.0/25) ;
  n14 = (1.0/25) ;
  n15 = (1.0/25) ;
  n16 = (1.0/25) ;
  n17 = (1.0/25) ;
  n18 = (1.0/25) ;
  n19 = (1.0/25) ;
  n20 = (1.0/25) ;
  n21 = (1.0/25) ;
  n22 = (1.0/25) ;
  n23 = (1.0/25) ;
  n24 = (1.0/25) ;
  
  // coordonnees absolues du point de base
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ; 
    
  outval0 = n0*tex2D(tex_img_inc, j-2, i-2 ) 
	+ n1*tex2D(tex_img_inc, j-1, i-2 ) 
	+ n2*tex2D(tex_img_inc, j  , i-2 ) 
	+ n3*tex2D(tex_img_inc, j+1, i-2 ) 
	+ n4*tex2D(tex_img_inc, j+2, i-2 ) 
    + n5*tex2D(tex_img_inc, j-2, i-1 ) 
	+ n6*tex2D(tex_img_inc, j-1, i-1 ) 
	+ n7*tex2D(tex_img_inc, j  , i-1 ) 
	+ n8*tex2D(tex_img_inc, j+1, i-1 ) 
	+ n9*tex2D(tex_img_inc, j+2, i-1 ) 
  + n10*tex2D(tex_img_inc, j-2, i ) 
  + n11*tex2D(tex_img_inc, j-1, i ) 
  + n12*tex2D(tex_img_inc, j  , i ) 
  + n13*tex2D(tex_img_inc, j+1, i ) 
  + n14*tex2D(tex_img_inc, j+2, i ) 
	+ n15*tex2D(tex_img_inc, j-2, i+1 ) 
	+ n16*tex2D(tex_img_inc, j-1, i+1 ) 
	+ n17*tex2D(tex_img_inc, j  , i+1 ) 
	+ n18*tex2D(tex_img_inc, j+1, i+1 ) 
	+ n19*tex2D(tex_img_inc, j+2, i+1 ) 
	+ n20*tex2D(tex_img_inc, j-2, i+2 ) 
  + n21*tex2D(tex_img_inc, j-1, i+2 ) 
  + n22*tex2D(tex_img_inc, j  , i+2 ) 
  + n23*tex2D(tex_img_inc, j+1, i+2 ) 
	+ n24*tex2D(tex_img_inc, j+2, i+2 ) ;
  
  output[ __mul24(i, j_dim) + j  ] = outval0 ;
}

// convolution non separable
// image en texture et noyau en registres
__global__ void kernel_convoGene7Reg8( unsigned char  *output, int i_dim, int j_dim)
{
  float outval0=0.0 ;
  float n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n13,n14,n15,n16,n17,n18,n19,n20,n21,n22,n23,n24 ;
  float n25,n26,n27,n28,n29,n30,n31,n32,n33,n34,n35,n36,n37,n38,n39,n40,n41,n42,n43,n44,n45,n46,n47,n48 ;
  
  n0 = (1.0/49) ;
  n1 = (1.0/49) ;
  n2 = (1.0/49) ;
  n3 = (1.0/49) ;
  n4 = (1.0/49) ;
  n5 = (1.0/49) ;
  n6 = (1.0/49) ;
  n7 = (1.0/49) ;
  n8 = (1.0/49) ;
  n9 = (1.0/49) ;
  n10 = (1.0/49) ;
  n11 = (1.0/49) ;
  n12 = (1.0/49) ;
  n13 = (1.0/49) ;
  n14 = (1.0/49) ;
  n15 = (1.0/49) ;
  n16 = (1.0/49) ;
  n17 = (1.0/49) ;
  n18 = (1.0/49) ;
  n19 = (1.0/49) ;
  n20 = (1.0/49) ;
  n21 = (1.0/49) ;
  n22 = (1.0/49) ;
  n23 = (1.0/49) ;
  n24 = (1.0/49) ;

  n25 = (1.0/49) ;
  n26 = (1.0/49) ;
  n27 = (1.0/49) ;
  n28 = (1.0/49) ;
  n29 = (1.0/49) ;
  n30 = (1.0/49) ;
  n31 = (1.0/49) ;
  n32 = (1.0/49) ;
  n33 = (1.0/49) ;
  n34 = (1.0/49) ;
  n35 = (1.0/49) ;
  n36 = (1.0/49) ;
  n37 = (1.0/49) ;
  n38 = (1.0/49) ;
  n39 = (1.0/49) ;
  n40 = (1.0/49) ;
  n41 = (1.0/49) ;
  n42 = (1.0/49) ;
  n43 = (1.0/49) ;
  n44 = (1.0/49) ;
  n45 = (1.0/49) ;
  n46 = (1.0/49) ;
  n47 = (1.0/49) ;
  n48 = (1.0/49) ;
  
  
  
  // coordonnees absolues du point de base
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ; 
    
  outval0 =
	  n0*tex2D(tex_img_inc, j-3, i-3 ) 
	+ n1*tex2D(tex_img_inc, j-2, i-3 ) 
	+ n2*tex2D(tex_img_inc, j-1, i-3 ) 
	+ n3*tex2D(tex_img_inc, j  , i-3 ) 
	+ n4*tex2D(tex_img_inc, j+1, i-3 ) 
    + n5*tex2D(tex_img_inc, j+2, i-3 ) 
	+ n6*tex2D(tex_img_inc, j+3, i-3 )
	+  n7*tex2D(tex_img_inc, j-3, i-2 ) 
	+  n8*tex2D(tex_img_inc, j-2, i-2 ) 
	+  n9*tex2D(tex_img_inc, j-1, i-2 ) 
	+ n10*tex2D(tex_img_inc, j  , i-2 ) 
	+ n11*tex2D(tex_img_inc, j+1, i-2 ) 
    + n12*tex2D(tex_img_inc, j+2, i-2 ) 
	+ n13*tex2D(tex_img_inc, j+3, i-2 )
	+ n14*tex2D(tex_img_inc, j-3, i-1 ) 
	+ n15*tex2D(tex_img_inc, j-2, i-1 ) 
	+ n16*tex2D(tex_img_inc, j-1, i-1 ) 
	+ n17*tex2D(tex_img_inc, j  , i-1 ) 
	+ n18*tex2D(tex_img_inc, j+1, i-1 ) 
    + n19*tex2D(tex_img_inc, j+2, i-1 ) 
	+ n20*tex2D(tex_img_inc, j+3, i-1 )
	+ n21*tex2D(tex_img_inc, j-3, i ) 
	+ n22*tex2D(tex_img_inc, j-2, i ) 
	+ n23*tex2D(tex_img_inc, j-1, i ) 
	+ n24*tex2D(tex_img_inc, j  , i ) 
	+ n25*tex2D(tex_img_inc, j+1, i ) 
    + n26*tex2D(tex_img_inc, j+2, i ) 
	+ n27*tex2D(tex_img_inc, j+3, i )
	+ n28*tex2D(tex_img_inc, j-3, i+1 ) 
	+ n29*tex2D(tex_img_inc, j-2, i+1 ) 
	+ n30*tex2D(tex_img_inc, j-1, i+1 ) 
	+ n31*tex2D(tex_img_inc, j  , i+1 ) 
	+ n32*tex2D(tex_img_inc, j+1, i+1 ) 
    + n33*tex2D(tex_img_inc, j+2, i+1 ) 
	+ n34*tex2D(tex_img_inc, j+3, i+1 )
	+ n35*tex2D(tex_img_inc, j-3, i+2 ) 
	+ n36*tex2D(tex_img_inc, j-2, i+2 ) 
	+ n37*tex2D(tex_img_inc, j-1, i+2 ) 
	+ n38*tex2D(tex_img_inc, j  , i+2 ) 
	+ n39*tex2D(tex_img_inc, j+1, i+2 ) 
    + n40*tex2D(tex_img_inc, j+2, i+2 ) 
	+ n41*tex2D(tex_img_inc, j+3, i+2 )
	+ n42*tex2D(tex_img_inc, j-3, i+3 ) 
	+ n43*tex2D(tex_img_inc, j-2, i+3 ) 
	+ n44*tex2D(tex_img_inc, j-1, i+3 ) 
	+ n45*tex2D(tex_img_inc, j  , i+3 ) 
	+ n46*tex2D(tex_img_inc, j+1, i+3 ) 
    + n47*tex2D(tex_img_inc, j+2, i+3 ) 
	+ n48*tex2D(tex_img_inc, j+3, i+3 ) ;
	
  
  output[ __mul24(i, j_dim) + j  ] = outval0 ;
}


__global__ void kernel_convoNonSep16( unsigned short*output, int i_dim, int j_dim)
{
  int idb, ic, jc ;
  int r = rnoyau ;
  int L = 2*r+1 ;
  int N = L*L ;
  float outval0=0 ;
    
  // coordonnees absolues du point de base
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ; 
  
  for (idb=0 ; idb< N ; idb++)
	{
	  ic = i-r + idb/L ;
	  jc = j-r + idb - (idb/L)*L ;
	  outval0 += ( noyau[ idb ]*tex2D(tex_img_in, jc, ic)) ;
	}

  output[ __mul24(i, j_dim) + j  ] = outval0 ;
}

// convolution non separable
// image en texture et noyau en mem constante
// prefetch des pixels en smem
// 1 pixel par thread 
__global__ void kernel_convoNonSepSh( unsigned short*output, int i_dim, int j_dim)
{
  int idb, ic, jc ;
  int L = 2*rnoyau+1 ;
  float outval0=0 ;
    
  // coordonnees absolues du point de base
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image

  // chargement en smem
  int idrow = threadIdx.y*(blockDim.x+2*rnoyau) ;
  
  extern __shared__ int roi[];

  // bloc 0 (en haut \E0 gauche) 
  roi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-rnoyau, i-rnoyau) ;
  // bloc 1 (en haut \E0 droite)...
  if ( threadIdx.x < L-1 ) //...ou plutot ce qu'il en manque
	roi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-rnoyau, i-rnoyau ) ;
  // bloc 2 ( en bas \E0 gauche)
  if ( threadIdx.y < L-1 )
	{
	  idrow = (threadIdx.y+blockDim.y)*(blockDim.x+L-1) ;
	  roi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-rnoyau, i+blockDim.y-rnoyau ) ;
	  //bloc 4 ( en bas \E0 doite )...
	  if ( threadIdx.x < L-1 ) //...ou ce qu'il en manque
		roi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-rnoyau, i+blockDim.y-rnoyau ) ;
	}
  __syncthreads();
  
  // calculs de convolution
  for (ic=0 ; ic<L ; ic++)
	for( jc=0 ; jc<L ; jc++)
   	  outval0 += ( noyau[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(blockDim.x+L-1)) + jc+threadIdx.x]) ;
	
  // 1 pixel par thread --> global mem
  output[ idx  ] = outval0 ;

}

__global__ void kernel_convoNonSepSh_2p(unsigned short *output, int i_dim, int j_dim)
{
  int idb, ic, jc ;
  int L = 2*rnoyau+1 ;
  float outval0=0.0, outval1=0.0 ;
  int bdimX = blockDim.x<<1 ;
  int tidX = threadIdx.x<<1 ;
    
  // coordonnees absolues du point de base
  int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image

  // chargement en smem
  int idrow = threadIdx.y*(bdimX+L-1) ;
  
  extern __shared__ int roi[];

  // bloc 0 (en haut \E0 gauche) 
  roi[ idrow  + tidX ] = tex2D(tex_img_ins, j-rnoyau, i-rnoyau) ;
  roi[ idrow  + tidX +1] = tex2D(tex_img_ins, j-rnoyau+1, i-rnoyau) ;

  // bloc 1 (en haut \E0 droite)...
  if ( threadIdx.x < rnoyau ) //...ou plutot ce qu'il en manque
	{
	  roi[ idrow + tidX+bdimX   ] = tex2D( tex_img_ins, j+bdimX-rnoyau, i-rnoyau ) ;
	  roi[ idrow + tidX+bdimX +1] = tex2D( tex_img_ins, j+bdimX-rnoyau+1, i-rnoyau ) ;
	}

  // bloc 2 ( en bas \E0 gauche)
  if ( threadIdx.y < L-1 )
	{
	  idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
	  roi[ idrow + tidX ] = tex2D( tex_img_ins, j-rnoyau, i+blockDim.y-rnoyau ) ;
	  roi[ idrow + tidX +1] = tex2D( tex_img_ins, j-rnoyau+1, i+blockDim.y-rnoyau ) ;
	  //bloc 4 ( en bas \E0 doite )...
	  
	  if ( 2*threadIdx.x < L-1 ) //...ou ce qu'il en manque
		{
		  roi[ idrow + tidX+bdimX ] = tex2D( tex_img_ins, j+bdimX-rnoyau, i+blockDim.y-rnoyau ) ;
		  roi[ idrow + tidX+bdimX +1] = tex2D( tex_img_ins, j+bdimX-rnoyau+1, i+blockDim.y-rnoyau ) ;
		  }
	}
  __syncthreads();
  
  // calculs de convolution
  for (ic=0 ; ic<L ; ic++)
	for( jc=0 ; jc<L ; jc++)
	  {
		outval0 += (noyau[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(2*blockDim.x+L-1)) + jc+2*threadIdx.x]) ;
		outval1 += (noyau[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(2*blockDim.x+L-1)) + jc+2*threadIdx.x+1]) ;																										 
	  }
	
  // 1 pixel par thread --> global mem
  output[ idx  ] = outval0 ;
  output[ idx+1 ] = outval1 ;
}

__global__ void kernel_convoNonSepSh_2p(unsigned char *output, int i_dim, int j_dim, int r)
{
  int idb, ic, jc ;
  int L = 2*r+1 ;
  float outval0=0.0, outval1=0.0 ;
  int bdimX = blockDim.x<<1 ;
  int tidX = threadIdx.x<<1 ;
    
  // coordonnees absolues du point de base
  int j = __umul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ; 
  int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image

  // chargement en smem
  int idrow = threadIdx.y*(bdimX+L-1) ;
  
  extern __shared__ int roi[];

  // bloc 0 (en haut \E0 gauche) 
  roi[ idrow  + tidX    ] = tex2D( tex_img_inc, j-r  , i-r) ;
  roi[ idrow  + tidX +1 ] = tex2D( tex_img_inc, j-r+1, i-r) ;

  // bloc 1 (en haut \E0 droite)...
  if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
	{
	  roi[ idrow + bdimX + tidX    ] = tex2D( tex_img_inc, j+bdimX-r  , i-r ) ;
	  roi[ idrow + bdimX + tidX +1 ] = tex2D( tex_img_inc, j+bdimX-r+1, i-r ) ;
	  }

  // bloc 2 ( en bas \E0 gauche)
  if ( threadIdx.y < L-1 )
	{
	  idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
	  roi[ idrow + tidX ] = tex2D( tex_img_inc, j-r, i+blockDim.y-r ) ;
	  roi[ idrow + tidX +1] = tex2D( tex_img_inc, j-r+1, i+blockDim.y-r ) ;
	  //bloc 4 ( en bas \E0 doite )...
	  
	  if ( threadIdx.x < r ) //...ou ce qu'il en manque
		{
		  roi[ idrow + tidX+bdimX   ] = tex2D( tex_img_inc, j+bdimX-r  , i+blockDim.y-r ) ;
		  roi[ idrow + tidX+bdimX +1] = tex2D( tex_img_inc, j+bdimX-r+1, i+blockDim.y-r ) ;
		  }
	}
  __syncthreads();
  
  // calculs de convolution
  for (ic=0 ; ic<L ; ic++)
	for( jc=0 ; jc<L ; jc++)
	  {
		outval0 += (masque[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX]) ;
		outval1 += (masque[ ic*L+jc ]*roi[ __mul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX+1]) ;																										 
	  }
	
  // 1 pixel par thread --> global mem
  output[ idx  ] = outval0 ;
  output[ idx+1 ] = outval1 ;
}

//4 pixels en ligne par thread
__global__ void kernel_convoNonSepSh_4p(unsigned char *output, int i_dim, int j_dim, int r)
{
  int idb, ic, jc ;
  int L = 2*r+1 ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
  int bdimX = blockDim.x<<2 ;
  int tidX = threadIdx.x<<2 ;
    
  // coordonnees absolues du point de base
  int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<2 ; 
  int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  int j0= __umul24(blockIdx.x,blockDim.x)<<2 ;
  
  int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image

  
  // chargement en smem
  int idrow = threadIdx.y*(bdimX+L-1) ;
  
  extern __shared__ unsigned char roi4p[];

  // bloc 0 (en haut \E0 gauche) 
  roi4p[ idrow  + tidX   ] = tex2D(tex_img_inc, j-r  , i-r) ;
  roi4p[ idrow  + tidX +1] = tex2D(tex_img_inc, j-r+1, i-r) ;
  roi4p[ idrow  + tidX +2] = tex2D(tex_img_inc, j-r+2, i-r) ;
  roi4p[ idrow  + tidX +3] = tex2D(tex_img_inc, j-r+3, i-r) ;
  
  // bloc 1 (en haut \E0 droite)...
  if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
	{
	  roi4p[ idrow + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r ) ;
	  roi4p[ idrow + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r ) ;
	}
  // bloc 2 ( en bas \E0 gauche)
  if ( threadIdx.y < L-1 )
	{
	  idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
	  roi4p[ idrow + tidX   ] = tex2D( tex_img_inc, j-r  , i+blockDim.y-r ) ;
	  roi4p[ idrow + tidX +1] = tex2D( tex_img_inc, j-r+1, i+blockDim.y-r ) ;
	  roi4p[ idrow + tidX +2] = tex2D( tex_img_inc, j-r+2, i+blockDim.y-r ) ;
	  roi4p[ idrow + tidX +3] = tex2D( tex_img_inc, j-r+3, i+blockDim.y-r ) ;
	  
	  //bloc 4 ( en bas \E0 doite )...
	  if ( threadIdx.x < r ) //...ou ce qu'il en manque
		{
		  roi4p[ idrow + bdimX +threadIdx.x   ] = tex2D( tex_img_inc, j0-r +bdimX +threadIdx.x, i+blockDim.y-r ) ;
		  roi4p[ idrow + bdimX +threadIdx.x +r] = tex2D( tex_img_inc, j0   +bdimX +threadIdx.x, i+blockDim.y-r ) ;
		}
	}
  __syncthreads();
  
  // calculs de convolution
  for (ic=0 ; ic<L ; ic++)
	for( jc=0 ; jc<L ; jc++)
	  {
		int baseRoi = __umul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX ;
		float valMask = masque[ __umul24(ic,L)+jc ] ;
		
		outval0 += valMask*roi4p[ baseRoi ] ;
		outval1 += valMask*roi4p[ baseRoi +1 ] ;
		outval2 += valMask*roi4p[ baseRoi +2 ] ;
		outval3 += valMask*roi4p[ baseRoi +3 ] ;
	  }
	
  // 1 pixel par thread --> global mem
  output[ idx   ] = outval0 ;
  output[ idx+1 ] = outval1 ;
  output[ idx+2 ] = outval2 ;
  output[ idx+3 ] = outval3 ;
}


//4 pixels en carre par thread
__global__ void kernel_convoNonSepSh_4pcarre(unsigned char *output, int i_dim, int j_dim, int r)
{
  int idb, ic, jc ;
  int L = 2*r+1 ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0 ;
  int bdimX = blockDim.x<<1 ;
  int tidX = threadIdx.x<<1 ;
  int bdimY = blockDim.y<<1 ;
  int tidY = threadIdx.y<<1 ;
    
  // coordonnees absolues du point de base
  int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<1 ; 
  int i = (__umul24( blockIdx.y, blockDim.y) + threadIdx.y)<<1 ;
  int j0= __umul24(blockIdx.x,blockDim.x)<<1 ;
  int i0= __umul24(blockIdx.y,blockDim.y)<<1 ;
  
  int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image

  
  // chargement en smem
  int idrowBase   = tidY*(bdimX+L-1) ;
  int idrowSecond = (tidY+1)*(bdimX+L-1) ;

  
  extern __shared__ unsigned char roi4p2[];

  // bloc 0 (en haut \E0 gauche) 
  roi4p2[ idrowBase  + tidX   ] = tex2D(tex_img_inc, j-r  , i-r  ) ;
  roi4p2[ idrowBase  + tidX +1] = tex2D(tex_img_inc, j-r+1, i-r  ) ;
  roi4p2[ idrowSecond+ tidX   ] = tex2D(tex_img_inc, j-r  , i-r+1) ;
  roi4p2[ idrowSecond+ tidX +1] = tex2D(tex_img_inc, j-r+1, i-r+1) ;
  
  // bloc 1 (en haut \E0 droite)...
  if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
	{
	  roi4p2[ idrowBase + bdimX   + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r ) ;
	  roi4p2[ idrowBase + bdimX   + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r ) ;
	  roi4p2[ idrowSecond + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r+1 ) ;
	  roi4p2[ idrowSecond + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r+1 ) ;
	}
  // bloc 2 ( en bas \E0 gauche)
  if ( threadIdx.y < L-1 )
	{
	  idrowBase   = (bdimY + threadIdx.y)*(bdimX+L-1) ;
	  roi4p2[ idrowBase + tidX     ] = tex2D( tex_img_inc, j-r  , i0-r +bdimY +threadIdx.y ) ;
	  roi4p2[ idrowBase + tidX   +1] = tex2D( tex_img_inc, j-r+1, i0-r +bdimY +threadIdx.y ) ;
	  	  
	  //bloc 4 ( en bas \E0 doite )...
	  if ( threadIdx.x < r ) //...ou ce qu'il en manque
		{
		  roi4p2[ idrowBase + bdimX + threadIdx.x    ]   = tex2D( tex_img_inc, j0-r +bdimX +threadIdx.x  , i0-r +bdimY +threadIdx.y ) ;
		  roi4p2[ idrowBase + bdimX + threadIdx.x +r ]   = tex2D( tex_img_inc, j0   +bdimX +threadIdx.x  , i0-r +bdimY +threadIdx.y ) ;
		}
	}
  __syncthreads();
  
  // calculs de convolution
  for (ic=0 ; ic<L ; ic++)
	for( jc=0 ; jc<L ; jc++)
	  {
		int baseRoi = __umul24(ic+tidY,(bdimX+L-1)) + jc+tidX ;
		float valMask = masque[ __umul24(ic,L)+jc ] ;
		
		outval0 += valMask*roi4p2[ baseRoi ] ;
		outval1 += valMask*roi4p2[ baseRoi +1 ] ;
		outval2 += valMask*roi4p2[ baseRoi +bdimX+L-1 ] ;
		outval3 += valMask*roi4p2[ baseRoi +bdimX+L ] ;
	  }
	
  // 1 pixel par thread --> global mem
  output[ idx   ] = outval0 ;
  output[ idx+1 ] = outval1 ;
  output[ idx+j_dim ] = outval2 ;
  output[ idx+j_dim+1 ] = outval3 ;
}


//8 pixels en ligne par thread
__global__ void kernel_convoNonSepSh_8p(unsigned char *output, int i_dim, int j_dim, int r)
{
  int idb, ic, jc;
  int L = 2*r+1 ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  int bdimX = blockDim.x<<3 ;
  int tidX = threadIdx.x<<3 ;
    
  // coordonnees absolues du point de base
  int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<3 ; 
  int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  int j0= __umul24(blockIdx.x,blockDim.x)<<3 ;
  
  int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image

  
  // chargement en smem
  int idrow = threadIdx.y*(bdimX+L-1) ;
  
  extern __shared__ unsigned char roi8p[];

  // bloc 0 (en haut \E0 gauche)
  for (int p=0; p<8; p++)
  roi8p[ idrow  + tidX +p ] = tex2D(tex_img_inc, j-r+p  , i-r) ;
  
  // bloc 1 (en haut \E0 droite)...
  if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
	{
	  roi8p[ idrow + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r ) ;
	  roi8p[ idrow + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r ) ;
	}
  // bloc 2 ( en bas \E0 gauche)
  if ( threadIdx.y < L-1 )
	{
	  idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
	  for (int p=0; p<8; p++)
		roi8p[ idrow + tidX +p  ] = tex2D( tex_img_inc, j-r+p  , i+blockDim.y-r ) ;
	  	  
	  //bloc 4 ( en bas \E0 doite )...
	  if ( threadIdx.x < r ) //...ou ce qu'il en manque
		{
		  roi8p[ idrow + bdimX +threadIdx.x   ] = tex2D( tex_img_inc, j0-r +bdimX +threadIdx.x, i+blockDim.y-r ) ;
		  roi8p[ idrow + bdimX +threadIdx.x +r] = tex2D( tex_img_inc, j0   +bdimX +threadIdx.x, i+blockDim.y-r ) ;
		}
	}
  __syncthreads();
  
  // calculs de convolution
  for (ic=0 ; ic<L ; ic++)
	for( jc=0 ; jc<L ; jc++)
	  {
		int baseRoi = __umul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX ;
		float valMask = masque[ __umul24(ic,L)+jc ] ;
		outval0 += valMask*roi8p[ baseRoi ] ;
		outval1 += valMask*roi8p[ baseRoi +1 ] ;
		outval2 += valMask*roi8p[ baseRoi +2 ] ;
		outval3 += valMask*roi8p[ baseRoi +3 ] ;
		outval4 += valMask*roi8p[ baseRoi +4 ] ;
		outval5 += valMask*roi8p[ baseRoi +5 ] ;
		outval6 += valMask*roi8p[ baseRoi +6 ] ;
		outval7 += valMask*roi8p[ baseRoi +7 ] ;
	  }
	
  // 1 pixel par thread --> global mem
  output[ idx   ] = outval0 ;
  output[ idx+1 ] = outval1 ;
  output[ idx+2 ] = outval2 ;
  output[ idx+3 ] = outval3 ;
  output[ idx+4 ] = outval4 ;
  output[ idx+5 ] = outval5 ;
  output[ idx+6 ] = outval6 ;
  output[ idx+7 ] = outval7 ;
}


//16 pixels en ligne par thread
__global__ void kernel_convoNonSepSh_16p(unsigned char *output, int i_dim, int j_dim, int r)
{
  int idb, ic, jc;
  int L = 2*r+1 ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  float outval8=0.0, outval9=0.0, outval10=0.0, outval11=0.0, outval12=0.0, outval13=0.0, outval14=0.0, outval15=0.0 ;
  int bdimX = blockDim.x<<4 ;
  int tidX = threadIdx.x<<4 ;
    
  // coordonnees absolues du point de base
  int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<4 ; 
  int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  int j0= __umul24(blockIdx.x,blockDim.x)<<4 ;
  
  int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image

  
  // chargement en smem
  int idrow = threadIdx.y*(bdimX+L-1) ;
  
  extern __shared__ unsigned char roi16p[];

  // bloc 0 (en haut \E0 gauche)
  for (int p=0; p<16; p++)
	roi16p[ idrow  + tidX +p ] = tex2D(tex_img_inc, j-r+p  , i-r) ;
  
  // bloc 1 (en haut \E0 droite)...
  if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
	{
	  roi16p[ idrow + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i-r ) ;
	  roi16p[ idrow + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i-r ) ;
	}
  // bloc 2 ( en bas \E0 gauche)
  if ( threadIdx.y < L-1 )
	{
	  idrow = (threadIdx.y+blockDim.y)*(bdimX+L-1) ;
	  for (int p=0; p<16; p++)
		roi16p[ idrow + tidX +p  ] = tex2D( tex_img_inc, j-r+p  , i+blockDim.y-r ) ;
	  	  
	  //bloc 4 ( en bas \E0 doite )...
	  if ( threadIdx.x < r ) //...ou ce qu'il en manque
		{
		  roi16p[ idrow + bdimX +threadIdx.x   ] = tex2D( tex_img_inc, j0-r +bdimX +threadIdx.x, i+blockDim.y-r ) ;
		  roi16p[ idrow + bdimX +threadIdx.x +r] = tex2D( tex_img_inc, j0   +bdimX +threadIdx.x, i+blockDim.y-r ) ;
		}
	}
  __syncthreads();
  
  // calculs de convolution
  for (ic=0 ; ic<L ; ic++)
	for( jc=0 ; jc<L ; jc++)
	  {
		int baseRoi = __umul24(ic+threadIdx.y,(bdimX+L-1)) + jc+tidX ;
		float valMask = masque[ __umul24(ic,L)+jc ] ;
		outval0 += valMask*roi16p[ baseRoi ] ;
		outval1 += valMask*roi16p[ baseRoi +1 ] ;
		outval2 += valMask*roi16p[ baseRoi +2 ] ;
		outval3 += valMask*roi16p[ baseRoi +3 ] ;
		outval4 += valMask*roi16p[ baseRoi +4 ] ;
		outval5 += valMask*roi16p[ baseRoi +5 ] ;
		outval6 += valMask*roi16p[ baseRoi +6 ] ;
		outval7 += valMask*roi16p[ baseRoi +7 ] ;
		outval8 += valMask*roi16p[ baseRoi +8] ;
		outval9 += valMask*roi16p[ baseRoi +9 ] ;
		outval10 += valMask*roi16p[ baseRoi +10 ] ;
		outval11 += valMask*roi16p[ baseRoi +11 ] ;
		outval12 += valMask*roi16p[ baseRoi +12 ] ;
		outval13 += valMask*roi16p[ baseRoi +13 ] ;
		outval14 += valMask*roi16p[ baseRoi +14 ] ;
		outval15 += valMask*roi16p[ baseRoi +15 ] ;
	  }
	
  // 1 pixel par thread --> global mem
  output[ idx   ] = outval0 ;
  output[ idx+1 ] = outval1 ;
  output[ idx+2 ] = outval2 ;
  output[ idx+3 ] = outval3 ;
  output[ idx+4 ] = outval4 ;
  output[ idx+5 ] = outval5 ;
  output[ idx+6 ] = outval6 ;
  output[ idx+7 ] = outval7 ;
  output[ idx+8 ] = outval8 ;
  output[ idx+9 ] = outval9 ;
  output[ idx+10 ] = outval10 ;
  output[ idx+11 ] = outval11 ;
  output[ idx+12 ] = outval12 ;
  output[ idx+13 ] = outval13 ;
  output[ idx+14 ] = outval14 ;
  output[ idx+15 ] = outval15 ;
}



// convolution non separable
// image en texture et noyau en mem constante
// fetch direct des pixels
// 2 pixels traités par thread => meilleur débit
__global__ void kernel_convoNonSep_2p( int *output, int i_dim, int j_dim)
{
  int idb, ic, jc ;
  int r = rnoyau ;
  int L = 2*r+1 ;
  int N = L*L ;
  float outval0=0, outval1=0 ;
    
  // coordonnees absolues du point de base
  int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x, 2) ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  int idx = __mul24(i, j_dim) + j ;              // indice  dans l'image

  #pragma unroll
  for (idb=0 ; idb< N ; idb++)
	{
	  ic = i-r + idb/L ;
	  jc = j-r + idb - (idb/L)*L ;
	  outval0 += ( noyau[ idb ]*tex2D(tex_img_in, jc, ic)) ;
	  outval1 += ( noyau[ idb ]*tex2D(tex_img_in, jc+1, ic)) ;
	}

  output[ idx   ] = outval0  ;
  output[ idx+1 ] = outval1  ;
  

}


// convolution non separable
// image en texture et noyau en mem constante
// fetch direct des pixels
// 2 pixels traités par thread => meilleur débit
__global__ void kernel_convoNonSep_2p_s( unsigned short*output, int i_dim, int j_dim)
{
  int idb, ic, jc ;
  int r = rnoyau ;
  int L = 2*r+1 ;
  int N = L*L ;
  float outval0=0, outval1=0 ;
      
  // coordonnees absolues du point de base
  int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x, 2) ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  int idx = __mul24(i, j_dim) + j ;              // indice  dans l'image

  #pragma unroll
  for (idb=0 ; idb< N ; idb++)
	{
	  ic = i-r + idb/L ;
	  jc = j-r + idb - (idb/L)*L ;
	  outval0 += ( noyau[ idb ]*tex2D(tex_img_ins, jc, ic)) ;
	  outval1 += ( noyau[ idb ]*tex2D(tex_img_ins, jc+1, ic)) ;
	}

  output[ idx   ] = outval0  ;
  output[ idx+1 ] = outval1  ;

}



// convolution separable
// image en texture et noyau en mem constante
__global__ void kernel_convoSep8V( unsigned char  *output, int i_dim, int j_dim, int r)
{
  int ic ;
  int L=2*r+1 ;
  float outval0=0.0 ;
  
  // coordonnees absolues du point de base
  int j = __mul24( blockIdx.x, blockDim.x ) + threadIdx.x ; 
  int i = __mul24( blockIdx.y, blockDim.y ) + threadIdx.y ; 

  //vertical
  for (ic=0 ; ic<L ; ic++)
	  outval0 += masque[ ic ]*tex2D(tex_img_inc, j, i-r+ic) ;

  output[ __mul24(i, j_dim) + j ] = outval0 ;
}




__global__ void kernel_convoSep8H( unsigned char  *output, int i_dim, int j_dim, int r)
{
  int jc ;
  int L=2*r+1 ;
  float outval0=0.0 ;
  
  // coordonnees absolues du point de base
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ; 

  //horizontal
  for (jc=0 ; jc<L ; jc++)
	  outval0 += masque[ L+jc ]*tex2D(tex_img_inc, j-r+jc, i) ;

  output[ __mul24(i, j_dim) + j ] = outval0 ;
}

// convolution separable
// image en texture et noyau en mem constante
__global__ void kernel_convoSep8Vx8p( unsigned char  *output, int i_dim, int j_dim, int r)
{
  int ic, baseIdx ;
  int L=2*r+1 ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  
  // coordonnees absolues du point de base
  int j = (__mul24( blockIdx.x, blockDim.x ) + threadIdx.x)<<3 ; 
  int i = __mul24( blockIdx.y, blockDim.y ) + threadIdx.y ; 

  //vertical
  for (ic=0 ; ic<L ; ic++)
	{
	  float valMask = masque[ ic ] ;
	  outval0 += valMask*tex2D(tex_img_inc, j  , i-r+ic) ;
	  outval1 += valMask*tex2D(tex_img_inc, j+1, i-r+ic) ;
	  outval2 += valMask*tex2D(tex_img_inc, j+2, i-r+ic) ;
	  outval3 += valMask*tex2D(tex_img_inc, j+3, i-r+ic) ;
	  outval4 += valMask*tex2D(tex_img_inc, j+4, i-r+ic) ;
	  outval5 += valMask*tex2D(tex_img_inc, j+5, i-r+ic) ;
	  outval6 += valMask*tex2D(tex_img_inc, j+6, i-r+ic) ;
	  outval7 += valMask*tex2D(tex_img_inc, j+7, i-r+ic) ;
	}

  baseIdx = __mul24(i  , j_dim) + j ;
  output[ baseIdx++ ] = outval0 ;
  output[ baseIdx++ ] = outval1 ;
  output[ baseIdx++ ] = outval2 ;
  output[ baseIdx++ ] = outval3 ;
  output[ baseIdx++ ] = outval4 ;
  output[ baseIdx++ ] = outval5 ;
  output[ baseIdx++ ] = outval6 ;
  output[ baseIdx++ ] = outval7 ;
}



// convolution separable
// image en texture et noyau en mem constante
__global__ void kernel_convoSep8HL3x8p( unsigned char  *output, int i_dim, int j_dim)
{
  int jc, baseIdx ;
  int L=3 ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  float val ;
  // coordonnees absolues du point de base
  int j = (__mul24( blockIdx.x, blockDim.x ) + threadIdx.x)<<3 ; 
  int i = __mul24( blockIdx.y, blockDim.y ) + threadIdx.y ; 

  //horizontal			
	val = tex2D(tex_img_inc, j-1  , i) ;
//p1 	
	outval0 += masque[0]*val ;
//p2
val = tex2D(tex_img_inc, j  , i) ;
	outval0 += masque[1]*val ;
	outval1 += masque[0]*val ;
//p3
val = tex2D(tex_img_inc, j+1  , i) ;
	outval0 += masque[2]*val ;
 	outval1 += masque[1]*val ;
        outval2 += masque[0]*val ;
//p4
val = tex2D(tex_img_inc, j+2  , i) ;
	outval1 += masque[2]*val ;
 	outval2 += masque[1]*val ;
        outval3 += masque[0]*val ;
//p5
val = tex2D(tex_img_inc, j+3  , i) ;
	outval2 += masque[2]*val ;
 	outval3 += masque[1]*val ;
        outval4 += masque[0]*val ;	
//p6
val = tex2D(tex_img_inc, j+4  , i) ;
	outval3 += masque[2]*val ;
 	outval4 += masque[1]*val ;
        outval5 += masque[0]*val ;
//p7
val = tex2D(tex_img_inc, j+5  , i) ;
	outval4 += masque[2]*val ;
 	outval5 += masque[1]*val ;
        outval6 += masque[0]*val ;
//p8
val = tex2D(tex_img_inc, j+6  , i) ;
	outval5 += masque[2]*val ;
 	outval6 += masque[1]*val ;
        outval7 += masque[0]*val ;
//p9
val = tex2D(tex_img_inc, j+7  , i) ;
	outval6 += masque[2]*val ;
 	outval7 += masque[1]*val ;
//p10
val = tex2D(tex_img_inc, j+8  , i) ;
	outval7 += masque[2]*val ;
 	
  baseIdx = __mul24(i  , j_dim) + j ;
  output[ baseIdx++ ] = outval0 ;
  output[ baseIdx++ ] = outval1 ;
  output[ baseIdx++ ] = outval2 ;
  output[ baseIdx++ ] = outval3 ;
  output[ baseIdx++ ] = outval4 ;
  output[ baseIdx++ ] = outval5 ;
  output[ baseIdx++ ] = outval6 ;
  output[ baseIdx++ ] = outval7 ;
}


// convolution separable
// image en texture et noyau en mem constante
__global__ void kernel_convoSep8HL3x8pG( unsigned char *input, unsigned char  *output, int i_dim, int j_dim)
{
  int jc, baseIdx, id ;
  int L=3 ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  float val ;
  // coordonnees absolues du point de base
  int j = (__mul24( blockIdx.x, blockDim.x ) + threadIdx.x)<<3 ; 
  int i = __mul24( blockIdx.y, blockDim.y ) + threadIdx.y ; 

  //horizontal	
baseIdx = __mul24(i  , j_dim) + j ;		
if (baseIdx >0 ) id = baseIdx-1 ; else id = baseIdx ;
	val = input[id] ;
//p1 	
	outval0 += masque[0]*val ;
//p2
val = input[baseIdx++] ;
	outval0 += masque[1]*val ;
	outval1 += masque[0]*val ;
//p3
val = input[baseIdx++] ;
	outval0 += masque[2]*val ;
 	outval1 += masque[1]*val ;
        outval2 += masque[0]*val ;
//p4
val = input[baseIdx++] ;
	outval1 += masque[2]*val ;
 	outval2 += masque[1]*val ;
        outval3 += masque[0]*val ;
//p5
val = input[baseIdx++] ;
	outval2 += masque[2]*val ;
 	outval3 += masque[1]*val ;
        outval4 += masque[0]*val ;	
//p6
val = input[baseIdx++] ;
	outval3 += masque[2]*val ;
 	outval4 += masque[1]*val ;
        outval5 += masque[0]*val ;
//p7
val = input[baseIdx++] ;
	outval4 += masque[2]*val ;
 	outval5 += masque[1]*val ;
        outval6 += masque[0]*val ;
//p8
val = input[baseIdx++] ;
	outval5 += masque[2]*val ;
 	outval6 += masque[1]*val ;
        outval7 += masque[0]*val ;
//p9
val = input[baseIdx++] ;
	outval6 += masque[2]*val ;
 	outval7 += masque[1]*val ;
//p10
val = input[baseIdx++] ;
	outval7 += masque[2]*val ;
 	
  baseIdx = __mul24(i  , j_dim) + j ;
  output[ baseIdx++ ] = outval0 ;
  output[ baseIdx++ ] = outval1 ;
  output[ baseIdx++ ] = outval2 ;
  output[ baseIdx++ ] = outval3 ;
  output[ baseIdx++ ] = outval4 ;
  output[ baseIdx++ ] = outval5 ;
  output[ baseIdx++ ] = outval6 ;
  output[ baseIdx++ ] = outval7 ;
}

// convolution separable
// image en texture et noyau en mem constante
__global__ void kernel_convoSep8HL5x8pG( unsigned char *input, unsigned char  *output, int i_dim, int j_dim)
{
  int jc, baseIdx, id ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  float val ;
  // coordonnees absolues du point de base
  int j = (__mul24( blockIdx.x, blockDim.x ) + threadIdx.x)<<3 ; 
  int i = __mul24( blockIdx.y, blockDim.y ) + threadIdx.y ; 

  //horizontal	
baseIdx = __mul24(i  , j_dim) + j ;		
if (baseIdx >1 ) id = baseIdx-2 ; else id = baseIdx ;
val = input[id] ;
//p1 	
	outval0 += masque[0]*val ;
//p2
if (baseIdx >0 ) id = baseIdx-1 ; else id = baseIdx ;
val = input[id] ;
	outval0 += masque[1]*val ;
	outval1 += masque[0]*val ;
//p3
val = input[baseIdx++] ;
	outval0 += masque[2]*val ;
 	outval1 += masque[1]*val ;
        outval2 += masque[0]*val ;
//p4
val = input[baseIdx++] ;
	outval0 += masque[3]*val ;
	outval1 += masque[2]*val ;
 	outval2 += masque[1]*val ;
        outval3 += masque[0]*val ;
//p5
val = input[baseIdx++] ;
	outval0 += masque[4]*val ;
	outval1 += masque[3]*val ;
	outval2 += masque[2]*val ;
 	outval3 += masque[1]*val ;
        outval4 += masque[0]*val ;	
//p6
val = input[baseIdx++] ;
	outval1 += masque[4]*val ;
	outval2 += masque[3]*val ;
	outval3 += masque[2]*val ;
 	outval4 += masque[1]*val ;
        outval5 += masque[0]*val ;
//p7
val = input[baseIdx++] ;
	outval2 += masque[4]*val ;
	outval3 += masque[3]*val ;
	outval4 += masque[2]*val ;
 	outval5 += masque[1]*val ;
        outval6 += masque[0]*val ;
//p8
val = input[baseIdx++] ;
	outval3 += masque[4]*val ;
	outval4 += masque[3]*val ;
	outval5 += masque[2]*val ;
 	outval6 += masque[1]*val ;
        outval7 += masque[0]*val ;
//p9
val = input[baseIdx++] ;
	outval4 += masque[4]*val ;
	outval5 += masque[3]*val ;
 	outval6 += masque[2]*val ;
        outval7 += masque[1]*val ;
//p10
val = input[baseIdx++] ;
	outval5 += masque[4]*val ;
 	outval6 += masque[3]*val ;
        outval7 += masque[2]*val ;
//p11
val = input[baseIdx++] ;
 	outval6 += masque[4]*val ;
        outval7 += masque[3]*val ;
//p11
val = input[baseIdx++] ;
        outval7 += masque[4]*val ;

 	
  baseIdx = __mul24(i  , j_dim) + j ;
  output[ baseIdx++ ] = outval0 ;
  output[ baseIdx++ ] = outval1 ;
  output[ baseIdx++ ] = outval2 ;
  output[ baseIdx++ ] = outval3 ;
  output[ baseIdx++ ] = outval4 ;
  output[ baseIdx++ ] = outval5 ;
  output[ baseIdx++ ] = outval6 ;
  output[ baseIdx++ ] = outval7 ;
}

// convolution separable
// image en texture et noyau en mem constante
__global__ void kernel_convoSep8HL7x8pG( unsigned char *input, unsigned char  *output, int i_dim, int j_dim)
{
  int jc, baseIdx, id ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  float val ;
  // coordonnees absolues du point de base
  int j = (__mul24( blockIdx.x, blockDim.x ) + threadIdx.x)<<3 ; 
  int i = __mul24( blockIdx.y, blockDim.y ) + threadIdx.y ; 

  //horizontal	
baseIdx = __mul24(i  , j_dim) + j ;		
if (baseIdx >2 ) id = baseIdx-3 ; else id = baseIdx ;
val = input[id] ;
//p1 	
	outval0 += masque[0]*val ;
//p2
if (baseIdx >1 ) id = baseIdx-2 ; else id = baseIdx ;
val = input[id] ;
	outval0 += masque[1]*val ;
	outval1 += masque[0]*val ;
//p3
if (baseIdx >0 ) id = baseIdx-1 ; else id = baseIdx ;
val = input[id] ;
	outval0 += masque[2]*val ;
 	outval1 += masque[1]*val ;
        outval2 += masque[0]*val ;
//p4
val = input[baseIdx++] ;
	outval0 += masque[3]*val ;
	outval1 += masque[2]*val ;
 	outval2 += masque[1]*val ;
        outval3 += masque[0]*val ;
//p5
val = input[baseIdx++] ;
	outval0 += masque[4]*val ;
	outval1 += masque[3]*val ;
	outval2 += masque[2]*val ;
 	outval3 += masque[1]*val ;
        outval4 += masque[0]*val ;	
//p6
val = input[baseIdx++] ;
	outval0 += masque[5]*val ;
	outval1 += masque[4]*val ;
	outval2 += masque[3]*val ;
 	outval3 += masque[2]*val ;
        outval4 += masque[1]*val ;
	outval5 += masque[0]*val ;
//p7
val = input[baseIdx++] ;
	outval0 += masque[6]*val ;
	outval1 += masque[5]*val ;
	outval2 += masque[4]*val ;
 	outval3 += masque[3]*val ;
        outval4 += masque[2]*val ;
	outval5 += masque[1]*val ;
	outval6 += masque[0]*val ;
//p8
val = input[baseIdx++] ;
	outval1 += masque[6]*val ;
	outval2 += masque[5]*val ;
	outval3 += masque[4]*val ;
 	outval4 += masque[3]*val ;
        outval5 += masque[2]*val ;
	outval6 += masque[1]*val ;
	outval7 += masque[0]*val ;
//p9
val = input[baseIdx++] ;
	outval2 += masque[6]*val ;
	outval3 += masque[5]*val ;
 	outval4 += masque[4]*val ;
        outval5 += masque[3]*val ;
	outval6 += masque[2]*val ;
	outval7 += masque[1]*val ;
//p10
val = input[baseIdx++] ;
	outval3 += masque[6]*val ;
 	outval4 += masque[5]*val ;
        outval5 += masque[4]*val ;
	outval6 += masque[3]*val ;
	outval7 += masque[2]*val ;
//p11
val = input[baseIdx++] ;
 	outval4 += masque[6]*val ;
        outval5 += masque[5]*val ;
	outval6 += masque[4]*val ;
	outval7 += masque[3]*val ;
//p12
val = input[baseIdx++] ;
        outval5 += masque[6]*val ;
	outval6 += masque[5]*val ;
	outval7 += masque[4]*val ;
//p13
val = input[baseIdx++] ;
       	outval6 += masque[6]*val ;
	outval7 += masque[5]*val ;
//p14
val = input[baseIdx++] ;
	outval7 += masque[6]*val ; 	

  baseIdx = __mul24(i  , j_dim) + j ;
  output[ baseIdx++ ] = outval0 ;
  output[ baseIdx++ ] = outval1 ;
  output[ baseIdx++ ] = outval2 ;
  output[ baseIdx++ ] = outval3 ;
  output[ baseIdx++ ] = outval4 ;
  output[ baseIdx++ ] = outval5 ;
  output[ baseIdx++ ] = outval6 ;
  output[ baseIdx++ ] = outval7 ;
}


//8 pixels en ligne par thread
__global__ void kernel_convoSepShx8pV(unsigned char *output, int i_dim, int j_dim, int r)
{
  int idb, ic, jc, p;
  int L = 2*r+1 ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  int bdimX = blockDim.x<<3 ;
  int tidX = threadIdx.x<<3 ;
    
  // coordonnees absolues du point de base
  int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<3 ; 
  int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  
  int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image

  
  // chargement en smem
  int idrow = threadIdx.y*bdimX ;
  
  extern __shared__ unsigned char roi8p[];

  // bloc 0 (en haut)
  for (p=0; p<8; p++)
	roi8p[ idrow  + tidX +p ] = tex2D(tex_img_inc, j+p  , i-r) ;
  
  
  // bloc 2 ( en bas)
  if ( threadIdx.y < L-1 )
	{
	  idrow = (threadIdx.y+blockDim.y)*bdimX ;
	  for (int p=0; p<8; p++)
		roi8p[ idrow + tidX +p  ] = tex2D( tex_img_inc, j+p  , i+blockDim.y-r ) ;
	  	  
	}
  __syncthreads();
  
  // calculs de convolution
  // passe verticale
  for (ic=0 ; ic<L ; ic++)
	  {
		int baseRoi = __umul24(ic+threadIdx.y,bdimX) + tidX ;
		float valMask = masque[ ic ] ;
		outval0 += valMask*roi8p[ baseRoi    ] ;
		outval1 += valMask*roi8p[ baseRoi +1 ] ;
		outval2 += valMask*roi8p[ baseRoi +2 ] ;
		outval3 += valMask*roi8p[ baseRoi +3 ] ;
		outval4 += valMask*roi8p[ baseRoi +4 ] ;
		outval5 += valMask*roi8p[ baseRoi +5 ] ;
		outval6 += valMask*roi8p[ baseRoi +6 ] ;
		outval7 += valMask*roi8p[ baseRoi +7 ] ;
	  }
	
  // 8 pixel par thread --> global mem
  output[ idx   ] = outval0 ;
  output[ idx+1 ] = outval1 ;
  output[ idx+2 ] = outval2 ;
  output[ idx+3 ] = outval3 ;
  output[ idx+4 ] = outval4 ;
  output[ idx+5 ] = outval5 ;
  output[ idx+6 ] = outval6 ;
  output[ idx+7 ] = outval7 ;
}

//8 pixels en ligne par thread
__global__ void kernel_convoSepShx8pVG(unsigned char *input, unsigned char *output, int i_dim, int j_dim, int r)
{
  int idb, ic, jc, p;
  int L = 2*r+1 ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  int bdimX = blockDim.x<<3 ;
  int tidX = threadIdx.x<<3 ;
    
  // coordonnees absolues du point de base
  int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<3 ; 
  int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  
  int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image

  
  // chargement en smem
  int idrow = threadIdx.y*bdimX ;
  
  extern __shared__ unsigned char roi8p[];

  // bloc 0 (en haut)
  for (p=0; p<8; p++)
	roi8p[ idrow  + tidX +p ] = input[ (i-r)*j_dim + j+p ] ;
  
  
  // bloc 2 ( en bas)
  if ( threadIdx.y < L-1 )
	{
	  idrow = (threadIdx.y+blockDim.y)*bdimX ;
	  for (int p=0; p<8; p++)
		roi8p[ idrow + tidX +p  ] = input[ (i+blockDim.y-r)*j_dim + j+p ] ;
	  	  
	}
  __syncthreads();
  
  // calculs de convolution
  // passe verticale
  for (ic=0 ; ic<L ; ic++)
	  {
		int baseRoi = __umul24(ic+threadIdx.y,bdimX) + tidX ;
		float valMask = masque[ ic ] ;
		outval0 += valMask*roi8p[ baseRoi    ] ;
		outval1 += valMask*roi8p[ baseRoi +1 ] ;
		outval2 += valMask*roi8p[ baseRoi +2 ] ;
		outval3 += valMask*roi8p[ baseRoi +3 ] ;
		outval4 += valMask*roi8p[ baseRoi +4 ] ;
		outval5 += valMask*roi8p[ baseRoi +5 ] ;
		outval6 += valMask*roi8p[ baseRoi +6 ] ;
		outval7 += valMask*roi8p[ baseRoi +7 ] ;
	  }
	
  // 8 pixel par thread --> global mem
  output[ idx   ] = outval0 ;
  output[ idx+1 ] = outval1 ;
  output[ idx+2 ] = outval2 ;
  output[ idx+3 ] = outval3 ;
  output[ idx+4 ] = outval4 ;
  output[ idx+5 ] = outval5 ;
  output[ idx+6 ] = outval6 ;
  output[ idx+7 ] = outval7 ;
}


__global__ void kernel_convoSepShx8pH(unsigned char *output, int i_dim, int j_dim, int r)
{
  int idb, ic, jc, p;
  int L = 2*r+1 ;
  float outval0=0.0, outval1=0.0, outval2=0.0, outval3=0.0, outval4=0.0, outval5=0.0, outval6=0.0, outval7=0.0 ;
  int bdimX = blockDim.x<<3 ;
  int tidX = threadIdx.x<<3 ;
    
  // coordonnees absolues du point de base
  int j = (__umul24(blockIdx.x,blockDim.x) + threadIdx.x)<<3 ; 
  int i = __umul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  int j0= __umul24(blockIdx.x,blockDim.x)<<3 ;  
  int idx = __umul24(i,j_dim) + j ;              // indice  dans l'image

  
  // chargement en smem
  int idrow = threadIdx.y*(bdimX+L-1) ;
  
  extern __shared__ unsigned char roi8p[];

  // bloc 0 (a gauche)
  for (p=0; p<8; p++)
	roi8p[  idrow  + tidX +p ] = tex2D(tex_img_inc, j-r+p  , i) ;

  // a droite
  if ( threadIdx.x < r ) //...ou plutot ce qu'il en manque
	{
	  roi8p[  idrow + bdimX + threadIdx.x    ] = tex2D( tex_img_inc, j0-r +bdimX+threadIdx.x  , i ) ;
	  roi8p[  idrow + bdimX + threadIdx.x +r ] = tex2D( tex_img_inc, j0   +bdimX+threadIdx.x  , i ) ;
	}
  
  __syncthreads();
  
  // calculs de convolution
  // passe horizontale
  for (jc=0 ; jc<L ; jc++)
	  {
		int baseRoi = idrow + tidX +jc ;
		float valMask = masque[ jc ] ;
		outval0 += valMask*roi8p[  baseRoi    ] ;
		outval1 += valMask*roi8p[  baseRoi +1 ] ;
		outval2 += valMask*roi8p[  baseRoi +2 ] ;
		outval3 += valMask*roi8p[  baseRoi +3 ] ;
		outval4 += valMask*roi8p[  baseRoi +4 ] ;
		outval5 += valMask*roi8p[  baseRoi +5 ] ;
		outval6 += valMask*roi8p[  baseRoi +6 ] ;
		outval7 += valMask*roi8p[  baseRoi +7 ] ;
	  }
	
  // 1 pixel par thread --> global mem
  output[ idx   ] = outval0 ;
  output[ idx+1 ] = outval1 ;
  output[ idx+2 ] = outval2 ;
  output[ idx+3 ] = outval3 ;
  output[ idx+4 ] = outval4 ;
  output[ idx+5 ] = outval5 ;
  output[ idx+6 ] = outval6 ;
  output[ idx+7 ] = outval7 ;
}


/*************************************************************************************************
 ***********************************************************************************************

              FIN DES kERNELS de CONVOLUTION

 ***********************************************************************************************
 *************************************************************************************************/
// kernel de la libjacket
// Exchange trick: Morgan McGuire, ShaderX 2008
#define s2(a,b)            { float tmp = a; a = min(a,b); b = max(tmp,b); }
#define mn3(a,b,c)         s2(a,b); s2(a,c);
#define mx3(a,b,c)         s2(b,c); s2(a,c);

#define mnmx3(a,b,c)       mx3(a,b,c); s2(a,b);                               // 3 exchanges
#define mnmx4(a,b,c,d)     s2(a,b); s2(c,d); s2(a,c); s2(b,d);                // 4 exchanges
#define mnmx5(a,b,c,d,e)   s2(a,b); s2(c,d); mn3(a,c,e); mx3(b,d,e);          // 6 exchanges
#define mnmx6(a,b,c,d,e,f) s2(a,d); s2(b,e); s2(c,f); mn3(a,b,c); mx3(d,e,f); // 7 exchanges

#define IN(X,Y)  ((0 <= (X) && (X) < nx && 0 <= (Y) && (Y) < ny) ? d_in[(Y)*nx+(X)] : 0)

__global__ static void kernel_medjacket(int nx, int ny, unsigned short*d_out, unsigned short*d_in)
{
  int x = __mul24(blockIdx.x , blockDim.x) + threadIdx.x;
  int y = __mul24(blockIdx.y , blockDim.y) + threadIdx.y;

    // pull top six from image
    float v[6] = { IN(x-1, y-1), IN(x  , y-1), IN(x+1, y-1),
                   IN(x-1, y  ), IN(x  , y  ), IN(x+1, y  ) };

    // with each pass, remove min and max values and add new value
    mnmx6(v[0], v[1], v[2], v[3], v[4], v[5]);
    v[5] = IN(x-1, y+1); // add new contestant
    mnmx5(v[1], v[2], v[3], v[4], v[5]);
    v[5] = IN(x  , y+1); // add new contestant
    mnmx4(v[2], v[3], v[4], v[5]);
    v[5] = IN(x+1, y+1); // add last contenstant
    mnmx3(v[3], v[4], v[5]);

    // pick the middle one
    d_out[y*nx + x] = v[4];
}

 
/***************************************************************
 *   fonction de tri de 2 valeurs entieres (min en premier)
 ***************************************************************/
__device__ inline void s(int* a, int* b)
{
  
  int tmp ;
  if (*a > *b)
	{
	  tmp = *b ;
	  *b = *a ;
	  *a = tmp ;
	  }
}

__device__ inline void s(unsigned short * a, unsigned short* b)
{
  
  unsigned short tmp ;
  if (*a > *b)
	{
	  tmp = *b ;
	  *b = *a ;
	  *a = tmp ;
	  }
}

__device__ inline void s(unsigned char * a, unsigned char * b)
{
  
  unsigned short tmp ;
  if (*a > *b)
	{
	  tmp = *b ;
	  *b = *a ;
	  *a = tmp ;
	  }
}

/***************************************************************
 *   fonction de min-max d'un tableau de n valeurs entieres 
 ***************************************************************/
__device__ void minmaxN(unsigned short* v, int n)
{
  for (int i=1; i< n; i++)
	s( v, v+i) ;
  for (int i=n-2; i>0; i--)
	s( v+i, v+n-1) ;
}

/***************************************************************
 *   fonction de tri naif d'un tableau de n valeurs entieres 
 ***************************************************************/
__device__ void bubTriN(int * v, int n)
{
  for (int i=0; i< n-1; i++)
	for (int j=i+1; j<n; j++)
	  s( v+i, v+j) ;
}


/*****************************************************************************************
 *    MACROS pour tri min-max utilisant la fonction de tri s(a,b) ci-dessus 
 *****************************************************************************************/
#define min3(a, b, c) s(a, b); s(a, c);
#define max3(a, b, c) s(b, c); s(a, c);
#define minmax3(a, b, c) max3(a, b, c); s(a, b); // max 3
#define minmax4(a, b, c, d) s(a, b); s(c, d); s(a, c); s(b, d); //max 4
#define minmax5(a, b, c, d, e) s(a, b); s(c, d); min3(a, c, e); max3(b, d, e); //max 6
#define minmax6(a, b, c, d, e, f) s(a,d); s(b, e); s(c, f); min3(a, b, c); max3(d, e, f); //max 7

/***************************************************************
 *   fonction de tri de 9 valeurs entieres 
 ***************************************************************/
#define bubTriReg9(a, b, c, d, e, f, g, h, i) s(a,b);s(c,d);s(e,f);minmax3(g,h,i);s(a,c);s(e,g);s(a,e);s(b,c);s(d,e);s(f,g);s(b,d);s(f,h);s(b,f);s(c,d);s(e,f);minmax3(g,h,i);s(c,e);s(c,g);s(d,e);s(f,g);s(d,f);s(d,h);s(e,f);minmax3(g,h,i);s(e,g);


__global__ void kernel_bubMedian3( unsigned short*output, int i_dim, int j_dim)
{  
  
  // coordonnees absolues du point
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
  
  /**************************************************************************
   *             tri(s)
   **************************************************************************/
  int a0, a1, a2, a3, a4, a5, a6, a7, a8 ;

  /********************************************************************************
   * affectation des valeurs
   ********************************************************************************/
  a0 = tex2D(tex_img_ins, j-1, i-1) ;
  a1 = tex2D(tex_img_ins, j  , i-1) ;
  a2 = tex2D(tex_img_ins, j+1, i-1) ;
  a3 = tex2D(tex_img_ins, j-1, i) ;
  a4 = tex2D(tex_img_ins, j  , i) ;
  a5 = tex2D(tex_img_ins, j+1, i) ;
  a6 = tex2D(tex_img_ins, j-1, i+1) ;
  a7 = tex2D(tex_img_ins, j  , i+1) ;
  a8 = tex2D(tex_img_ins, j+1, i+1) ;

  //tri selection
  bubTriReg9(&a0, &a1, &a2, &a3, &a4, &a5, &a6, &a7, &a8);
  
  //median au milieu !
  output[ __mul24(i, j_dim) +j ] = a4 ;
}

__global__ void kernel_median3( unsigned short*output, int i_dim, int j_dim)
{  
  
  // coordonnees absolues du point
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
  
  /**************************************************************************
   *             tri(s)
   **************************************************************************/
  int a0, a1, a2, a3, a4, a5 ;

  /********************************************************************************
   * les six premieres valeurs (suffisant pour median 3x3 par forgetfull selection)
   ********************************************************************************/
  a0 = tex2D(tex_img_ins, j-1, i-1) ;
  a1 = tex2D(tex_img_ins, j, i-1) ;
  a2 = tex2D(tex_img_ins, j+1, i-1) ;
  a3 = tex2D(tex_img_ins, j-1, i) ;
  a4 = tex2D(tex_img_ins, j, i) ;
  a5 = tex2D(tex_img_ins, j+1, i) ;

  
  //min max aux extremites
  minmax6(&a0, &a1, &a2, &a3, &a4, &a5) ;

  /********************************************
   * les deux valeurs suivantes aux extremites 
   ********************************************/
  a5 = tex2D(tex_img_in, j-1, i+1) ;

  minmax5(&a1, &a2, &a3, &a4, &a5) ;

  /********************************************
   * la derniere valeur a la fin
   ********************************************/
  
  a5 = tex2D(tex_img_ins, j, i+1) ;

  minmax4(&a2, &a3, &a4, &a5) ;

  a5 = tex2D(tex_img_ins, j+1, i+1) ;
  minmax3(&a3, &a4, &a5) ;
  
  
  //median au milieu !
  output[ __mul24(i, j_dim) +j ] = a4 ;
 
}

__global__ void kernel_median3Sh( unsigned short*output, int i_dim, int j_dim)
{  
  int idb, ic, jc ;
  int L = 3 ;
    
  // coordonnees absolues du point de base
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24( blockIdx.y, blockDim.y) + threadIdx.y ;
  int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image

  // chargement en smem
  int idrow = threadIdx.y*(blockDim.x+2) ;
    
  extern __shared__ int medRoi[];

  // bloc 0 (en haut \E0 gauche) 
  medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-1, i-1) ;
  // bloc 1 (en haut \E0 droite)...
  if ( threadIdx.x < L-1 ) //...ou plutot ce qu'il en manque
	medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-1, i-1 ) ;
  // bloc 2 ( en bas \E0 gauche)
  if ( threadIdx.y < L-1 )
	{
	  idrow = (threadIdx.y+blockDim.y)*(blockDim.x+L-1) ;
	  medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-1, i+blockDim.y-1 ) ;
	  //bloc 4 ( en bas \E0 doite )...
	  if ( threadIdx.x < L-1 ) //...ou ce qu'il en manque
		medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-1, i+blockDim.y-1 ) ;
	}
  __syncthreads();

  /**************************************************************************
   *             tri(s)
   **************************************************************************/
  int a0, a1, a2, a3, a4, a5 ;

  /********************************************************************************
   * les six premieres valeurs (suffisant pour median 3x3 par forgetfull selection)
   ********************************************************************************/
  // l'index dans medRoi est du type (threadIdx.y+ic)*(blockDim.X+L-1)+ (threadIdx.x+jc)
  // ou ic,jc sont les coordonnees du point dans le noyau 3x3
  a0 = medRoi[ (threadIdx.y)*(blockDim.x+L-1)+ (threadIdx.x+1) ] ; 
  a1 = medRoi[ (threadIdx.y)*(blockDim.x+L-1)+ (threadIdx.x+2)  ] ; 
  a2 = medRoi[ (threadIdx.y+1)*(blockDim.x+L-1)+ (threadIdx.x+1)] ; 
  a3 = medRoi[ (threadIdx.y+1)*(blockDim.x+L-1)+ (threadIdx.x+2)] ; 
  a4 = medRoi[ (threadIdx.y+2)*(blockDim.x+L-1)+ (threadIdx.x+1)] ; 
  a5 = medRoi[ (threadIdx.y+2)*(blockDim.x+L-1)+ (threadIdx.x+2)] ; 

  
  //min max aux extremites
  minmax6(&a0, &a1, &a2, &a3, &a4, &a5) ;

  /********************************************
   * les deux valeurs suivantes aux extremites 
   ********************************************/
  a5 = medRoi[ (threadIdx.y)*(blockDim.x+L-1)+ (threadIdx.x)] ; 

  minmax5(&a1, &a2, &a3, &a4, &a5) ;

  /********************************************
   * la derniere valeur a la fin
   ********************************************/
  
  a5 = medRoi[ (threadIdx.y+1)*(blockDim.x+L-1)+ (threadIdx.x)] ;

  minmax4(&a2, &a3, &a4, &a5) ;

  a5 = medRoi[ (threadIdx.y+2)*(blockDim.x+L-1)+ (threadIdx.x)] ;
  minmax3(&a3, &a4, &a5) ;
  
  
  //median au milieu !
  output[ __mul24(i, j_dim) +j ] = a4 ;
 
}

/************************************************
 * version qui gere 2 pxels par thread
 * Objectif plus déebit en cachant la latence
 * minimiser les acces en texture
 * defaut = plus de registres
 ************************************************/
__global__ void kernel_median3_2pix( unsigned short*output, int i_dim, int j_dim)
{  
  
  // coordonnees absolues du point
  int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
  
  /**************************************************************************
   *             tri(s)
   **************************************************************************/
  int a0, a1, a2, a3, a4, a5 ;
  int b0, b1, b2, b3, b4, b5 ;

  /********************************************************************************
   * les six premieres valeurs (suffisant pour median 3x3 par forgetfull selection)
   ********************************************************************************/
  a0 = tex2D(tex_img_ins, j  , i-1) ;
  a1 = tex2D(tex_img_ins, j+1, i-1) ;
  a2 = tex2D(tex_img_ins, j  , i  ) ;
  a3 = tex2D(tex_img_ins, j+1, i  ) ;
  a4 = tex2D(tex_img_ins, j  , i+1) ;
  a5 = tex2D(tex_img_ins, j+1, i+1) ;

  
  //min max aux extremites
  minmax6(&a0, &a1, &a2, &a3, &a4, &a5) ;

  b0=a0; b1=a1; b2=a2; b3=a3; b4=a4; b5=a5;
  
  /********************************************
   * la valeur suivante au bout 
   ********************************************/
  a5 = tex2D(tex_img_ins, j-1, i  ) ;
  b5 = tex2D(tex_img_ins, j+2, i  ) ;
  
  minmax5(&a1, &a2, &a3, &a4, &a5) ;
  minmax5(&b1, &b2, &b3, &b4, &b5) ;

  /********************************************
   * la derniere valeur a la fin
   ********************************************/
  a5 = tex2D(tex_img_ins, j-1, i-1) ;
  b5 = tex2D(tex_img_ins, j+2, i-1) ;
  
  minmax4(&a2, &a3, &a4, &a5) ;
  minmax4(&b2, &b3, &b4, &b5) ;
  
  a5 = tex2D(tex_img_ins, j-1, i+1) ;
  b5 = tex2D(tex_img_ins, j+2, i+1) ;
  
  minmax3(&a3, &a4, &a5) ;
  minmax3(&b3, &b4, &b5) ;
  
  //median au milieu !
  output[ __mul24(i, j_dim) +j ] = a4 ;
  output[ __mul24(i, j_dim) +j+1 ] = b4 ;
 
}

__global__ void kernel_median3_2pix8( unsigned char  *output, int i_dim, int j_dim)
{  
  
  // coordonnees absolues du point
  int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
  
  /**************************************************************************
   *             tri(s)
   **************************************************************************/
  int a0, a1, a2, a3, a4, a5 ;
  int b0, b1, b2, b3, b4, b5 ;

  /********************************************************************************
   * les six premieres valeurs (suffisant pour median 3x3 par forgetfull selection)
   ********************************************************************************/
  a0 = tex2D(tex_img_inc, j  , i-1) ;
  a1 = tex2D(tex_img_inc, j+1, i-1) ;
  a2 = tex2D(tex_img_inc, j  , i  ) ;
  a3 = tex2D(tex_img_inc, j+1, i  ) ;
  a4 = tex2D(tex_img_inc, j  , i+1) ;
  a5 = tex2D(tex_img_inc, j+1, i+1) ;

  
  //min max aux extremites
  minmax6(&a0, &a1, &a2, &a3, &a4, &a5) ;

  b0=a0; b1=a1; b2=a2; b3=a3; b4=a4; b5=a5;
  
  /********************************************
   * la valeur suivante au bout 
   ********************************************/
  a5 = tex2D(tex_img_inc, j-1, i  ) ;
  b5 = tex2D(tex_img_inc, j+2, i  ) ;
  
  minmax5(&a1, &a2, &a3, &a4, &a5) ;
  minmax5(&b1, &b2, &b3, &b4, &b5) ;

  /********************************************
   * la derniere valeur a la fin
   ********************************************/
  a5 = tex2D(tex_img_inc, j-1, i-1) ;
  b5 = tex2D(tex_img_inc, j+2, i-1) ;
  
  minmax4(&a2, &a3, &a4, &a5) ;
  minmax4(&b2, &b3, &b4, &b5) ;
  
  a5 = tex2D(tex_img_inc, j-1, i+1) ;
  b5 = tex2D(tex_img_inc, j+2, i+1) ;
  
  minmax3(&a3, &a4, &a5) ;
  minmax3(&b3, &b4, &b5) ;
  
  //median au milieu !
  output[ __mul24(i, j_dim) +j ] = a4 ;
  output[ __mul24(i, j_dim) +j+1 ] = b4 ;
 
}


#define minmax7(a,b,c,d,e,f,g) s(a,b);s(c,d);s(e,f);s(a,c);s(e,g);s(a,e);s(b,d);s(f,g);s(d,g);//max 9
#define minmax8(a,b,c,d,e,f,g,h) s(a,b);s(c,d);s(e,f);s(g,h);s(a,c);s(e,g);s(a,e);s(b,d);s(f,h);s(d,h);//max 10
#define minmax9(a,b,c,d,e,f,g,h,i) s(a,b);s(c,d);s(e,f);minmax3(g,h,i);s(a,c);s(e,g);s(a,e);s(b,d);s(f,i);s(d,i);//max 12
#define minmax10(a,b,c,d,e,f,g,h,i,j) s(a,b);s(c,d);s(e,f);s(g,h);s(i,j);s(a,c);min3(e,g,i);s(a,e);s(b,d);max3(f,h,j);s(d,j);//max 13
#define minmax11(a,b,c,d,e,f,g,h,i,j,k) s(a,b);s(c,d);s(e,f);s(g,h);min3(i,j,k);s(a,c);min3(e,g,i);s(a,e);s(b,d);max3(f,h,k);s(d,k);//max 14
#define minmax12(a,b,c,d,e,f,g,h,i,j,k,l) s(a,b);s(c,d);s(e,f);s(g,h);s(i,j);s(k,l);s(a,c);s(e,g);s(i,k);min3(a,e,i);s(b,d);s(f,h);s(j,l);max3(d,h,l); //max 16
#define minmax13(a,b,c,d,e,f,g,h,i,j,k,l,m) s(a,b);s(c,d);s(e,f);s(g,h);s(i,j);min3(k,l,m);s(a,c);s(e,g);s(i,k);min3(a,e,i);s(b,d);s(f,h);max3(j,l,m);max3(d,h,m); // max 18 !!
#define min14(a,b,c,d,e,f,g,h,i,j,k,l,m,n) s(a,b);s(c,d);s(e,f);s(g,h);s(i,j);s(k,l);s(m,n);s(a,c);s(e,g);s(i,k);s(a,e);s(i,m);s(a,i); // 13 mouvts max
#define max14(a,b,c,d,e,f,g,h,i,j,k,l,m,n) s(a,b);s(c,d);s(e,f);s(g,h);s(i,j);s(k,l);s(m,n);s(b,d);s(f,h);s(j,l);s(d,h);s(l,n);s(h,n); // 13 mouvts max
#define minmax14(a,b,c,d,e,f,g,h,i,j,k,l,m,n) min14(a,b,c,d,e,f,g,h,i,j,k,l,m,n); s(b,d);s(f,h);s(j,l);s(d,h);s(l,n);s(h,n);           // 19 mouvts max 

#define m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o) s(a,b);s(c,d);s(e,f);s(g,h);s(i,j);s(k,l);s(a,c);s(e,g);s(i,k);s(b,d);s(f,h);s(j,l); // max 12
#define minmax15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);minmax3(m,n,o);s(a,e);s(i,m);s(d,h);s(l,o);s(a,i);s(h,o); //max 21 
#define minmax16(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(m,o);s(a,e);s(i,m);s(n,p);s(d,h);s(l,p);s(a,i);s(h,p);//max 22
#define minmax17(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);minmax3(o,p,q);s(a,e);s(d,h);s(i,m);min3(a,i,o);max3(h,n,q);//max 24
#define minmax18(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);minmax3(m,o,q);minmax3(n,p,r);s(a,e);s(m,n);s(q,r);s(i,m);s(d,h);s(l,r);s(a,i);s(h,r); //max 29 !!!
#define minmax19(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);minmax3(q,r,ss);s(m,o);s(n,p);s(a,e);s(i,m);s(d,h);s(l,p);min3(a,i,q);max3(h,p,ss);//max 28
#define minmax20(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss,t) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);s(ss,t);s(m,o);s(q,ss);s(n,p);s(r,t);s(a,e);s(i,m);s(d,h);s(l,p);min3(a,i,q);max3(h,p,t);//max 29
#define minmax21(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss,t,u) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);minmax3(ss,t,u);s(m,o);s(q,ss);s(n,p);s(r,u);s(a,e);s(i,m);s(d,h);s(l,p);min3(a,i,q);max3(h,p,u); //max 30
#define minmax22(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss,t,u,v) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);s(ss,t);s(u,v);s(m,o);min3(q,ss,u);s(n,p);max3(r,t,v);s(a,e);s(i,m);s(d,h);s(l,p);min3(a,i,q);max3(h,p,v); //max31
#define minmax23(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss,t,u,v,w) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);s(ss,t);minmax3(u,v,w);s(m,o);s(q,ss);s(n,p);s(r,t);s(a,e);s(i,m);s(q,u);s(d,h);s(l,p);s(t,w);min3(a,i,q);max3(h,p,w); //max 33
#define minmax24(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss,t,u,v,w,x) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);s(ss,t);s(u,v);s(w,x);s(m,o);s(q,ss);s(u,w);s(n,p);s(r,t);s(v,x);s(a,e);s(i,m);s(q,u);s(d,h);s(l,p);s(t,x);min3(a,i,q);max3(h,p,x); //max 34
#define minmax25(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss,t,u,v,w,x,y) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);s(ss,t);s(u,v);minmax3(w,x,y);s(m,o);s(q,ss);s(u,w);s(n,p);s(r,t);s(v,y);s(a,e);s(i,m);s(q,u);s(d,h);s(l,p);s(t,y);min3(a,i,q);max3(h,p,y); //max 36
#define minmax26(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,ss,t,u,v,w,x,y,z) m15(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o);s(m,n);s(o,p);s(q,r);s(ss,t);s(u,v);s(w,x);s(y,z);s(m,o);s(q,ss);min3(u,w,y);s(n,p);s(r,t);max3(v,x,z);s(a,e);s(i,m);s(q,u);s(d,h);s(l,p);s(t,z);min3(a,i,q);max3(h,p,z);//max 37


__global__ void kernel_median5( unsigned short*output, int i_dim, int j_dim)
{  
  
  // coordonnees absolues du point
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
  
  /**************************************************************************
   *             tri(s)
   **************************************************************************/
  int a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13 ;

  /********************************************************************************
   * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
   ********************************************************************************/
  //premiere ligne
  a0  = tex2D(tex_img_ins, j-2, i-2) ;
  a1  = tex2D(tex_img_ins, j-1, i-2) ;
  a2  = tex2D(tex_img_ins, j  , i-2) ;
  a3  = tex2D(tex_img_ins, j+1, i-2) ;
  a4  = tex2D(tex_img_ins, j+2, i-2) ;
  //deuxieme ligne
  a5  = tex2D(tex_img_ins, j-2, i-1) ;
  a6  = tex2D(tex_img_ins, j-1, i-1) ;
  a7  = tex2D(tex_img_ins, j  , i-1) ;
  a8  = tex2D(tex_img_ins, j+1, i-1) ;
  a9  = tex2D(tex_img_ins, j+2, i-1) ;
  //les 4 premiers de la troisieme ligne
  a10 = tex2D(tex_img_ins, j-2, i) ;
  a11 = tex2D(tex_img_ins, j-1, i) ;
  a12 = tex2D(tex_img_ins, j  , i) ;
  a13 = tex2D(tex_img_ins, j+1, i) ;

  //min max aux extremites
  minmax14(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeurs suivante (15)
  a13 = tex2D(tex_img_ins, j+2, i);
  //minmax aux extremites
  minmax13(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  
  //chargement valeur suivante (16)
  a13 = tex2D(tex_img_ins, j-2, i+1);
  //minmax aux extremites
  minmax12(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (17)
  a13 = tex2D(tex_img_ins, j-1, i+1);
  //minmax aux extremites
  minmax11(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (18)
  a13 = tex2D(tex_img_ins, j  , i+1);
  //minmax aux extremites
  minmax10(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (19)
  a13 = tex2D(tex_img_ins, j+1, i+1);
  //minmax aux extremites
  minmax9(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (20)
  a13 = tex2D(tex_img_ins, j+2, i+1);
  //minmax aux extremites
  minmax8(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (21)
  a13 = tex2D(tex_img_ins, j-2, i+2);
  //minmax aux extremites
  minmax7(&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (22)
  a13 = tex2D(tex_img_ins, j-1, i+2);
  //minmax aux extremites
  minmax6(&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (23)
  a13 = tex2D(tex_img_ins, j  , i+2);
  //minmax aux extremites
  minmax5(&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (24)
  a13 = tex2D(tex_img_ins, j+1  , i+2);
  //minmax aux extremites
  minmax4(&a10,&a11,&a12,&a13);

  //chargement valeur suivante (25)
  a13 = tex2D(tex_img_ins, j+2  , i+2);
  //minmax aux extremites
  minmax3(&a11,&a12,&a13);

  
  //median au milieu !
  output[ __mul24(i, j_dim) +j ] = a12 ;
 
}

__global__ void kernel_median5_sh( unsigned short*output, int i_dim, int j_dim)
{  
  
  // coordonnees absolues du point
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
  int L = 5 ;
  int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image

  // chargement en smem
  int idrow = threadIdx.y*(blockDim.x+4) ;
    
  extern __shared__ int medRoi[];

  // bloc 0 (en haut \E0 gauche) 
  medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-2, i-2) ;
  // bloc 1 (en haut \E0 droite)...
  if ( threadIdx.x < L-1 ) //...ou plutot ce qu'il en manque
	medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-2, i-2 ) ;
  // bloc 2 ( en bas \E0 gauche)
  if ( threadIdx.y < L-1 )
	{
	  idrow = (threadIdx.y+blockDim.y)*(blockDim.x+L-1) ;
	  medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-2, i+blockDim.y-2 ) ;
	  //bloc 4 ( en bas \E0 doite )...
	  if ( threadIdx.x < L-1 ) //...ou ce qu'il en manque
		medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-2, i+blockDim.y-2 ) ;
	}
  __syncthreads();

  /**************************************************************************
   *             tri(s)
   **************************************************************************/
  int a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13 ;

  /********************************************************************************
   * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
   ********************************************************************************/
  // l'index dans medRoi est du type (threadIdx.y+ic)*(blockDim.x+L-1)+ (threadIdx.x+jc)
  // ou ic,jc sont les coordonnees du point dans le noya
  //premiere ligne
  a0  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i-2) ;
  a1  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i-2) ;
  a2  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i-2) ;
  a3  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i-2) ;
  a4  = medRoi[ (threadIdx.y )*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i-2) ;
  //deuxieme ligne
  a5  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i-1) ;
  a6  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i-1) ;
  a7  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i-1) ;
  a8  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i-1) ;
  a9  = medRoi[ (threadIdx.y +1)*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i-1) ;
  //les 4 premiers de la troisieme ligne
  a10 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i) ;
  a11 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i) ;
  a12 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i) ;
  a13 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i) ;

  //min max aux extremites
  minmax14(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeurs suivante (15)
  a13 = medRoi[ (threadIdx.y +2)*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i);
  //minmax aux extremites
  minmax13(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  
  //chargement valeur suivante (16)
  a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i+1);
  //minmax aux extremites
  minmax12(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (17)
  a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i+1);
  //minmax aux extremites
  minmax11(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (18)
  a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i+1);
  //minmax aux extremites
  minmax10(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (19)
  a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i+1);
  //minmax aux extremites
  minmax9(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (20)
  a13 = medRoi[ (threadIdx.y +3)*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i+1);
  //minmax aux extremites
  minmax8(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (21)
  a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i+2);
  //minmax aux extremites
  minmax7(&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (22)
  a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i+2);
  //minmax aux extremites
  minmax6(&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (23)
  a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i+2);
  //minmax aux extremites
  minmax5(&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (24)
  a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1  , i+2);
  //minmax aux extremites
  minmax4(&a10,&a11,&a12,&a13);

  //chargement valeur suivante (25)
  a13 = medRoi[ (threadIdx.y +4)*(blockDim.x+L-1)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2  , i+2);
  //minmax aux extremites
  minmax3(&a11,&a12,&a13);

  
  //median au milieu !
  output[ __mul24(i, j_dim) +j ] = a12 ;
 
}


__global__ void kernel_median5_sh_loc( unsigned short*output, int i_dim, int j_dim)
{  
  
  // coordonnees absolues du point
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
  // int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image

  // chargement en smem
  int idrow = threadIdx.y*(blockDim.x+4) ;
    
  extern __shared__ int medRoi[];

  // bloc 0 (en haut \E0 gauche) 
  medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-2, i-2) ;
  // bloc 1 (en haut \E0 droite)...
  if ( threadIdx.x < 4 ) //...ou plutot ce qu'il en manque
	medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-2, i-2 ) ;
  // bloc 2 ( en bas \E0 gauche)
  if ( threadIdx.y < 4 )
	{
	  idrow = (threadIdx.y+blockDim.y)*(blockDim.x+4) ;
	  medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-2, i+blockDim.y-2 ) ;
	  //bloc 4 ( en bas \E0 doite )...
	  if ( threadIdx.x < 4 ) //...ou ce qu'il en manque
		medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-2, i+blockDim.y-2 ) ;
	}
  __syncthreads();

  /**************************************************************************
   *             tri(s)
   **************************************************************************/
  /********************************************************************************
   * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
   ********************************************************************************/
  // l'index dans medRoi est du type (threadIdx.y+ic)*(blockDim.x+L-1)+ (threadIdx.x+jc)
  // ou ic,jc sont les coordonnees du point dans le noyau
  
  int a[] = { medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x ) ],
			  medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x +1) ],
			  medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x +2) ],
			  medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x +3) ],
			  medRoi[ (threadIdx.y )*(blockDim.x+4)+ (threadIdx.x +4) ],
			  medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x ) ],
			  medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x +1) ],
			  medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x +2) ],
			  medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x +3) ],
			  medRoi[ (threadIdx.y +1)*(blockDim.x+4)+ (threadIdx.x +4) ],
			  medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x ) ],
			  medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x +1) ],
			  medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x +2) ],
			  medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x +3) ]};


  //min max aux extremites
  minmax14(&a[0],&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  //chargement valeurs suivante (15)
  a[13] = medRoi[ (threadIdx.y +2)*(blockDim.x+4)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i);
  //minmax aux extremites
  minmax13(&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  
  //chargement valeur suivante (16)
  a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i+1);
  //minmax aux extremites
  minmax12(&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  //chargement valeur suivante (17)
  a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i+1);
  //minmax aux extremites
  minmax11(&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  //chargement valeur suivante (18)
  a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i+1);
  //minmax aux extremites
  minmax10(&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  //chargement valeur suivante (19)
  a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1, i+1);
  //minmax aux extremites
  minmax9(&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  //chargement valeur suivante (20)
  a[13] = medRoi[ (threadIdx.y +3)*(blockDim.x+4)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2, i+1);
  //minmax aux extremites
  minmax8(&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  //chargement valeur suivante (21)
  a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x ) ] ; //tex2D(tex_img_ins, j-2, i+2);
  //minmax aux extremites
  minmax7(&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  //chargement valeur suivante (22)
  a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x +1) ] ; //tex2D(tex_img_ins, j-1, i+2);
  //minmax aux extremites
  minmax6(&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  //chargement valeur suivante (23)
  a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x +2) ] ; //tex2D(tex_img_ins, j  , i+2);
  //minmax aux extremites
  minmax5(&a[9],&a[10],&a[11],&a[12],&a[13]);

  //chargement valeur suivante (24)
  a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x +3) ] ; //tex2D(tex_img_ins, j+1  , i+2);
  //minmax aux extremites
  minmax4(&a[10],&a[11],&a[12],&a[13]);

  //chargement valeur suivante (25)
  a[13] = medRoi[ (threadIdx.y +4)*(blockDim.x+4)+ (threadIdx.x +4) ] ; //tex2D(tex_img_ins, j+2  , i+2);
  //minmax aux extremites
  minmax3(&a[11],&a[12],&a[13]);

  
  //median au milieu !
  output[ __mul24(i, j_dim) +j ] = a[12] ;
 
}

__global__ void kernel_median5_sh_loc_2p( unsigned short*output, int i_dim, int j_dim)
{  
  
  // coordonnees absolues du point
  int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
  // int idx = __mul24(i,j_dim) + j ;              // indice  dans l'image
  int bdimX = blockDim.x<<1 ;
  int tidX = threadIdx.x<<1 ;
  
  // chargement en smem
  int idrow = threadIdx.y*(bdimX+2*2) ;
    
  extern __shared__ int medRoi[];

  // bloc 0 (en haut \E0 gauche) 
  medRoi[ idrow  + tidX ] = tex2D(tex_img_ins, j-2, i-2) ;
  medRoi[ idrow  + tidX +1] = tex2D(tex_img_ins, j-1, i-2) ;
  // bloc 1 (en haut \E0 droite)...
  if ( tidX < 4 ) //...ou plutot ce qu'il en manque
	{
	  medRoi[ idrow + tidX + bdimX ] = tex2D( tex_img_ins, j+bdimX-2, i-2 ) ;
	  medRoi[ idrow + tidX + bdimX +1] = tex2D( tex_img_ins, j+bdimX-1, i-2 ) ;
	}
  // bloc 2 ( en bas \E0 gauche)
  if ( threadIdx.y < 4 )
	{
	  idrow = (threadIdx.y+blockDim.y)*(bdimX+4) ;
	  medRoi[ idrow + tidX ] = tex2D( tex_img_ins, j-2, i+blockDim.y-2 ) ;
	  medRoi[ idrow + tidX +1] = tex2D( tex_img_ins, j-1, i+blockDim.y-2 ) ;
	  //bloc 4 ( en bas \E0 doite )...
	  if ( tidX < 4 ) //...ou ce qu'il en manque
		{
		  medRoi[ idrow + tidX + bdimX ] = tex2D( tex_img_ins, j+bdimX-2, i+blockDim.y-2 ) ;
		  medRoi[ idrow + tidX + bdimX +1] = tex2D( tex_img_ins, j+bdimX-1, i+blockDim.y-2 ) ;
		}
	}
  __syncthreads();

  /**************************************************************************
   *             tri(s)
   **************************************************************************/
  /********************************************************************************
   * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
   ********************************************************************************/
  // l'index dans medRoi est du type (threadIdx.y+ic)*(bdimX+L-1)+ (tidX+jc)
  // ou ic,jc sont les coordonnees du point dans le noyau

  //on commence le tri avec les pixels du voisinage commun aux 2 pixels \E0 calculer
  
  int a[] = { medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +1) ],
			  medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +2) ],
			  medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +3) ],
			  medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +4) ],
			  
			  medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +1) ],
			  medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +2) ],
			  medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +3) ],
			  medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +4) ],
			  
			  medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +1) ],
			  medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +2) ],
			  medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +3) ],
			  medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +4) ],
			  
			  medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +1) ],
			  medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +2) ]
  };

  //min max aux extremites
  minmax14(&a[0],&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  //chargement valeurs suivante (15)
  a[13] = medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +3) ] ; 
  //minmax aux extremites
  minmax13(&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  
  //chargement valeur suivante (16)
  a[13] = medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +4) ] ; 
  //minmax aux extremites
  minmax12(&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  //chargement valeur suivante (17)
  a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +1) ] ; 
  //minmax aux extremites
  minmax11(&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  //chargement valeur suivante (18)
  a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +2) ] ; 
  //minmax aux extremites
  minmax10(&a[4],&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  //chargement valeur suivante (19)
  a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +3) ] ; 
  //minmax aux extremites
  minmax9(&a[5],&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  //chargement valeur suivante (20)
  a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +4) ] ; 
  //minmax aux extremites
  minmax8(&a[6],&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);

  // fin des pixels communs. la suite doit se faire en 2 tris distincts.
  int b[] = { a[7], a[8], a[9], a[10], a[11], a[12], a[13]};

  //chargement valeur suivante (21)
  a[13] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX ) ] ;
  b[6 ] = medRoi[ (threadIdx.y +4)*(bdimX+4)+ (tidX +5) ] ;
  //minmax aux extremites
  minmax7(&a[7],&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
  minmax7(&b[0],&b[1],&b[2],&b[3],&b[4],&b[5],&b[6]);

  //chargement valeur suivante (22)
  a[13] = medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX ) ] ;
  b[6 ] = medRoi[ (threadIdx.y +3)*(bdimX+4)+ (tidX +5) ] ; 
  //minmax aux extremites
  minmax6(&a[8],&a[9],&a[10],&a[11],&a[12],&a[13]);
  minmax6(&b[1],&b[2],&b[3],&b[4],&b[5],&b[6]);

  //chargement valeur suivante (23)
  a[13] = medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX ) ] ;
  b[6 ] = medRoi[ (threadIdx.y +2)*(bdimX+4)+ (tidX +5) ] ; 
  //minmax aux extremites
  minmax5(&a[9],&a[10],&a[11],&a[12],&a[13]);
  minmax5(&b[2],&b[3],&b[4],&b[5],&b[6]);

  //chargement valeur suivante (24)
  a[13] = medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX ) ] ;
  b[6 ] = medRoi[ (threadIdx.y +1)*(bdimX+4)+ (tidX +5) ] ;
  
  //minmax aux extremites
  minmax4(&a[10],&a[11],&a[12],&a[13]);
  minmax4(&b[3],&b[4],&b[5],&b[6]);

  //chargement valeur suivante (25)
  a[13] = medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX ) ] ;
  b[6 ] = medRoi[ (threadIdx.y )*(bdimX+4)+ (tidX +5) ] ; 
  //minmax aux extremites
  minmax3(&a[11],&a[12],&a[13]);
  minmax3(&b[4],&b[5],&b[6]);
  
  //median au milieu !
  output[ __mul24(i, j_dim) +j ] = a[12] ;
  output[ __mul24(i, j_dim) +j+1 ] = b[5] ;
 
}

 
__global__ void kernel_medianR( unsigned short * output, int i_dim, int j_dim, int r)
{   
  // coordonnees absolues du point
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ; 
  unsigned short cpt ;
  unsigned short ic, jc ;
  
  unsigned short hist[256] ;
  for (ic =0; ic<256; ic++) hist[ic]=0; // init histogramme
  //for (ic =0; ic<25; ic++) hist[ic]=1; // init histogramme

  //generation histogramme
  for(ic=i-r; ic<=i+r; ic++ )
	for(jc=j-r; jc<=j+r; jc++)
	  hist[tex2D(tex_img_ins, jc, ic)]++ ;

  //parcours histogramme 
  cpt = 0 ;
  for(ic=0; ic<256; ic++)
	{
	  cpt += hist[ic] ;
	  //selection de la valeur du percentile (ici 50%=SUM/2)
	  if ( cpt > ((2*r+1)*(2*r+1))>>1 ) break ; 
	}
  
  output[ __mul24(i, j_dim) +j ] =  ic ;  
}

__global__ void kernel_medianRSH( unsigned short * output, int i_dim, int j_dim, int r)
{   
  // coordonnees absolues du point
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ; 
  unsigned short cpt ;
  unsigned short ic, jc ;

  // chargement en smem
  int idrow = threadIdx.y*(blockDim.x+2*r) ;
  
  extern __shared__ int medRoi[];

  // bloc 0 (en haut \E0 gauche) 
  medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-r, i-r) ;
  // bloc 1 (en haut \E0 droite)...
  if ( threadIdx.x < 2*r ) //...ou plutot ce qu'il en manque
	medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-r, i-r ) ;
  // bloc 2 ( en bas \E0 gauche)
  if ( threadIdx.y < 2*r )
	{
	  idrow = (threadIdx.y+blockDim.y)*(blockDim.x+2*r) ;
	  medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-r, i+blockDim.y-r ) ;
	  //bloc 4 ( en bas \E0 doite )...
	  if ( threadIdx.x < 2*r ) //...ou ce qu'il en manque
		medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-r, i+blockDim.y-r ) ;
	}
  __syncthreads();
  
  int hist[256] ;
  for (ic =0; ic<256; ic++) hist[ic]=0; // init histogramme

  //generation histogramme
  for(ic=0; ic<=2*r; ic++ )
	for(jc=0; jc<=2*r; jc++)
	  hist[medRoi[(threadIdx.y+ic)*(blockDim.x+2*r)+ (threadIdx.x+jc)]]++ ;

  //parcours histogramme 
  cpt = 0 ;
  for(ic=0; ic<256; ic++)
	{
	  cpt += hist[ic] ;
	  //selection de la valeur du percentile (ici 50%=SUM/2)
	  if ( cpt > (((2*r+1)*(2*r+1))>>1) ) break ; 
	}
  
  output[ __mul24(i, j_dim) +j ] =  ic ;  
}


__global__ void kernel_median5_2pix( unsigned short*output, int i_dim, int j_dim)
{  
  
  // coordonnees absolues du point
  int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
  
  /**************************************************************************
   *             tri(s)
   **************************************************************************/
  int a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13 ;
  int b7, b8, b9, b10, b11, b12, b13 ;															

  /********************************************************************************
   * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
   ********************************************************************************/
  //premiere ligne
  a0  = tex2D(tex_img_ins, j-1, i-2) ;
  a1  = tex2D(tex_img_ins, j  , i-2) ;
  a2  = tex2D(tex_img_ins, j+1, i-2) ;
  a3  = tex2D(tex_img_ins, j+2, i-2) ;
  //deuxieme ligne
  a4  = tex2D(tex_img_ins, j-1, i-1) ;
  a5  = tex2D(tex_img_ins, j  , i-1) ;
  a6  = tex2D(tex_img_ins, j+1, i-1) ;
  a7  = tex2D(tex_img_ins, j+2, i-1) ;
  //troisieme ligne
  a8  = tex2D(tex_img_ins, j-1, i) ;
  a9  = tex2D(tex_img_ins, j  , i) ;
  a10 = tex2D(tex_img_ins, j+1, i) ;
  a11 = tex2D(tex_img_ins, j+2, i) ;
  //les 2 premiers de la 4eme ligne
  a12 = tex2D(tex_img_ins, j-1, i+1) ;
  a13 = tex2D(tex_img_ins, j  , i+1) ;

  //min max aux extremites
  minmax14(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
  
  //chargement valeurs suivante (15)
  a13 = tex2D(tex_img_ins, j+1, i+1);
  //minmax aux extremites
  minmax13(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
  
  //chargement valeur suivante (16)
  a13 = tex2D(tex_img_ins, j+2, i+1);
  //minmax aux extremites
  minmax12(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (17)
  a13 = tex2D(tex_img_ins, j-1, i+2);
  //minmax aux extremites
  minmax11(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (18)
  a13 = tex2D(tex_img_ins, j  , i+2);
  //minmax aux extremites
  minmax10(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (19)
  a13 = tex2D(tex_img_ins, j+1, i+2);
  //minmax aux extremites
  minmax9(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (20)
  a13 = tex2D(tex_img_ins, j+2, i+2);
  //minmax aux extmites
  minmax8(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
  
  // fin des pixels voisins communs deux pixels centraux
  b7=a7; b8=a8; b9=a9; b10=a10; b11=a11; b12=a12; b13=a13;
  
  //chargement valeur suivante (21)
  a13 = tex2D(tex_img_ins, j-2, i-2);
  b13 = tex2D(tex_img_ins, j+3, i-2);
  //minmax aux extremites
  minmax7(&a7,&a8,&a9,&a10,&a11,&a12,&a13);
  minmax7(&b7,&b8,&b9,&b10,&b11,&b12,&b13);
  
  //chargement valeur suivante (22)
  a13 = tex2D(tex_img_ins, j-2, i-1);
  b13 = tex2D(tex_img_ins, j+3, i-1);
  //minmax aux extremites
  minmax6(&a8,&a9,&a10,&a11,&a12,&a13);
  minmax6(&b8,&b9,&b10,&b11,&b12,&b13);

  //chargement valeur suivante (23)
  a13 = tex2D(tex_img_ins, j-2, i  );
  b13 = tex2D(tex_img_ins, j+3, i  );
  //minmax aux extremites
  minmax5(&a9,&a10,&a11,&a12,&a13);
  minmax5(&b9,&b10,&b11,&b12,&b13);

  //chargement valeur suivante (24)
  a13 = tex2D(tex_img_ins, j-2, i+1);
  b13 = tex2D(tex_img_ins, j+3, i+1);
  //minmax aux extremites
  minmax4(&a10,&a11,&a12,&a13);
  minmax4(&b10,&b11,&b12,&b13);
  
  //chargement valeur suivante (25)
  a13 = tex2D(tex_img_ins, j-2, i+2);
  b13 = tex2D(tex_img_ins, j+3, i+2);
  //minmax aux extremites
  minmax3(&a11,&a12,&a13);
  minmax3(&b11,&b12,&b13);

  //median au milieu !
  output[ __mul24(i, j_dim) +j ] = a12 ;
  output[ __mul24(i, j_dim) +j+1 ] = b12 ;
 
}


__global__ void kernel_median5_2pix( unsigned char  *output, int i_dim, int j_dim)
{  
  
  // coordonnees absolues du point
  int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
  
  /**************************************************************************
   *             tri(s)
   **************************************************************************/
  int a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13 ;
  int b7, b8, b9, b10, b11, b12, b13 ;															

  /********************************************************************************
   * les 14 premieres valeurs (suffisant pour median 5x5 par forgetfull selection)
   ********************************************************************************/
  //premiere ligne
  a0  = tex2D(tex_img_inc, j-1, i-2) ;
  a1  = tex2D(tex_img_inc, j  , i-2) ;
  a2  = tex2D(tex_img_inc, j+1, i-2) ;
  a3  = tex2D(tex_img_inc, j+2, i-2) ;
  //deuxieme ligne
  a4  = tex2D(tex_img_inc, j-1, i-1) ;
  a5  = tex2D(tex_img_inc, j  , i-1) ;
  a6  = tex2D(tex_img_inc, j+1, i-1) ;
  a7  = tex2D(tex_img_inc, j+2, i-1) ;
  //troisieme ligne
  a8  = tex2D(tex_img_inc, j-1, i) ;
  a9  = tex2D(tex_img_inc, j  , i) ;
  a10 = tex2D(tex_img_inc, j+1, i) ;
  a11 = tex2D(tex_img_inc, j+2, i) ;
  //les 2 premiers de la 4eme ligne
  a12 = tex2D(tex_img_inc, j-1, i+1) ;
  a13 = tex2D(tex_img_inc, j  , i+1) ;

  //min max aux extremites
  minmax14(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
  
  //chargement valeurs suivante (15)
  a13 = tex2D(tex_img_inc, j+1, i+1);
  //minmax aux extremites
  minmax13(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
  
  //chargement valeur suivante (16)
  a13 = tex2D(tex_img_inc, j+2, i+1);
  //minmax aux extremites
  minmax12(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (17)
  a13 = tex2D(tex_img_inc, j-1, i+2);
  //minmax aux extremites
  minmax11(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (18)
  a13 = tex2D(tex_img_inc, j  , i+2);
  //minmax aux extremites
  minmax10(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (19)
  a13 = tex2D(tex_img_inc, j+1, i+2);
  //minmax aux extremites
  minmax9(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);

  //chargement valeur suivante (20)
  a13 = tex2D(tex_img_inc, j+2, i+2);
  //minmax aux extmites
  minmax8(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13);
  
  // fin des pixels voisins communs deux pixels centraux
  b7=a7; b8=a8; b9=a9; b10=a10; b11=a11; b12=a12; b13=a13;
  
  //chargement valeur suivante (21)
  a13 = tex2D(tex_img_inc, j-2, i-2);
  b13 = tex2D(tex_img_inc, j+3, i-2);
  //minmax aux extremites
  minmax7(&a7,&a8,&a9,&a10,&a11,&a12,&a13);
  minmax7(&b7,&b8,&b9,&b10,&b11,&b12,&b13);
  
  //chargement valeur suivante (22)
  a13 = tex2D(tex_img_inc, j-2, i-1);
  b13 = tex2D(tex_img_inc, j+3, i-1);
  //minmax aux extremites
  minmax6(&a8,&a9,&a10,&a11,&a12,&a13);
  minmax6(&b8,&b9,&b10,&b11,&b12,&b13);

  //chargement valeur suivante (23)
  a13 = tex2D(tex_img_inc, j-2, i  );
  b13 = tex2D(tex_img_inc, j+3, i  );
  //minmax aux extremites
  minmax5(&a9,&a10,&a11,&a12,&a13);
  minmax5(&b9,&b10,&b11,&b12,&b13);

  //chargement valeur suivante (24)
  a13 = tex2D(tex_img_inc, j-2, i+1);
  b13 = tex2D(tex_img_inc, j+3, i+1);
  //minmax aux extremites
  minmax4(&a10,&a11,&a12,&a13);
  minmax4(&b10,&b11,&b12,&b13);
  
  //chargement valeur suivante (25)
  a13 = tex2D(tex_img_inc, j-2, i+2);
  b13 = tex2D(tex_img_inc, j+3, i+2);
  //minmax aux extremites
  minmax3(&a11,&a12,&a13);
  minmax3(&b11,&b12,&b13);

  //median au milieu !
  output[ __mul24(i, j_dim) +j ] = a12 ;
  output[ __mul24(i, j_dim) +j+1 ] = b12 ;
 
}


__global__ void kernel_median7( unsigned short*output, int i_dim, int j_dim)
{  
  
  // coordonnees absolues du point
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
  
  /**************************************************************************
   *             tri(s)
   **************************************************************************/
  int a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,a25 ;												

  /****************************************
   * les 26 premieres valeurs 
   ****************************************/
  //premiere ligne
  a0  = tex2D(tex_img_ins, j-2, i-3) ;
  a1  = tex2D(tex_img_ins, j-1, i-3) ;
  a2  = tex2D(tex_img_ins, j  , i-3) ;
  a3  = tex2D(tex_img_ins, j+1, i-3) ;
  a4  = tex2D(tex_img_ins, j+2, i-3) ;
  a5  = tex2D(tex_img_ins, j+3, i-3) ;
  //deuxieme ligne
  a6  = tex2D(tex_img_ins, j-2, i-2) ;
  a7  = tex2D(tex_img_ins, j-1, i-2) ;
  a8  = tex2D(tex_img_ins, j  , i-2) ;
  a9  = tex2D(tex_img_ins, j+1, i-2) ;
  a10 = tex2D(tex_img_ins, j+2, i-2) ;
  a11 = tex2D(tex_img_ins, j+3, i-2) ;
  //troisieme ligne
  a12 = tex2D(tex_img_ins, j-2, i-1) ;
  a13 = tex2D(tex_img_ins, j-1, i-1) ;
  a14 = tex2D(tex_img_ins, j  , i-1) ;
  a15 = tex2D(tex_img_ins, j+1, i-1) ;
  a16 = tex2D(tex_img_ins, j+2, i-1) ;
  a17 = tex2D(tex_img_ins, j+3, i-1) ;
  //quatrieme ligne
  a18 = tex2D(tex_img_ins, j-2, i  ) ;
  a19 = tex2D(tex_img_ins, j-1, i  ) ;
  a20 = tex2D(tex_img_ins, j  , i  ) ;
  a21 = tex2D(tex_img_ins, j+1, i  ) ;
  a22 = tex2D(tex_img_ins, j+2, i  ) ;
  a23 = tex2D(tex_img_ins, j+3, i  ) ;
  //cinquieme ligne
  a24 = tex2D(tex_img_ins, j-2, i+1) ;
  a25 = tex2D(tex_img_ins, j-1, i+1) ;
  
  //min max aux extremites
  minmax26(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeurs suivante (26)
  a25 = tex2D(tex_img_ins, j , i+1);
  //minmax aux extremites
  minmax25(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (27)
  a25 = tex2D(tex_img_ins, j+1, i+1);
  //minmax aux extremites
  minmax24(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (28)
  a25 = tex2D(tex_img_ins, j+2, i+1);
  //minmax aux extremites
  minmax23(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (29)
  a25 = tex2D(tex_img_ins, j+3, i+1);
  //minmax aux extremites
  minmax22(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (30)
  a25 = tex2D(tex_img_ins, j-2, i+2);
  //minmax aux extremites
  minmax21(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (31)
  a25 = tex2D(tex_img_ins, j-1, i+2);
  //minmax aux extmites
  minmax20(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (32)
  a25 = tex2D(tex_img_ins, j  , i+2);
  //minmax aux extmites
  minmax19(&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (33)
  a25 = tex2D(tex_img_ins, j+1, i+2);
  //minmax aux extmites
  minmax18(&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (34)
  a25 = tex2D(tex_img_ins, j+2, i+2);
  //minmax aux extmites
  minmax17(&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (35)
  a25 = tex2D(tex_img_ins, j+3, i+2);
  //minmax aux extmites
  minmax16(&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (36)
  a25 = tex2D(tex_img_ins, j-2, i+3);
  //minmax aux extmites
  minmax15(&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (37)
  a25 = tex2D(tex_img_ins, j-1, i+3);
  //minmax aux extmites
  
  minmax14(&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  //chargement valeur suivante (38)
  a25 = tex2D(tex_img_ins, j  , i+3);
  //minmax aux extmites
  minmax13(&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (39)
  a25 = tex2D(tex_img_ins, j+1, i+3);
  //minmax aux extmites
  minmax12(&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (40)
  a25 = tex2D(tex_img_ins, j+2, i+3);
  //minmax aux extmites
  minmax11(&a15, &a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (41)
  a25 = tex2D(tex_img_ins, j+3, i+3);
  //minmax aux extmites
  minmax10(&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
    
  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_ins, j-3, i+3);
  //minmax aux extremites
  minmax9(&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_ins, j-3, i+2);
  //minmax aux extremites
  minmax8(&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_ins, j-3, i+1);
  //minmax aux extremites
  minmax7(&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_ins, j-3, i);
  //minmax aux extremites
  minmax6(&a20,&a21,&a22,&a23,&a24,&a25);
 
  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_ins, j-3, i-1);
  //minmax aux extremites
  minmax5(&a21,&a22,&a23,&a24,&a25);
 
  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_ins, j-3, i-2);
  //minmax aux extremites
  minmax4(&a22,&a23,&a24,&a25);
 
  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_ins, j-3, i-3);
  //minmax aux extremites
  minmax3(&a23,&a24,&a25);
  
  //medians au milieu !
  output[ __mul24(i, j_dim) +j ] = a24 ;
 
}



__global__ void kernel_median7_2pix( unsigned short*output, int i_dim, int j_dim)
{  
  
  // coordonnees absolues du point
  int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
  
  /**************************************************************************
   *             tri(s)
   **************************************************************************/
  int a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,a25 ;
  int b17,b18,b19,b20,b21,b22,b23,b24,b25;															

  /****************************************
   * les 26 premieres valeurs 
   ****************************************/
  //premiere ligne
  a0  = tex2D(tex_img_ins, j-2, i-3) ;
  a1  = tex2D(tex_img_ins, j-1, i-3) ;
  a2  = tex2D(tex_img_ins, j  , i-3) ;
  a3  = tex2D(tex_img_ins, j+1, i-3) ;
  a4  = tex2D(tex_img_ins, j+2, i-3) ;
  a5  = tex2D(tex_img_ins, j+3, i-3) ;
  //deuxieme ligne
  a6  = tex2D(tex_img_ins, j-2, i-2) ;
  a7  = tex2D(tex_img_ins, j-1, i-2) ;
  a8  = tex2D(tex_img_ins, j  , i-2) ;
  a9  = tex2D(tex_img_ins, j+1, i-2) ;
  a10 = tex2D(tex_img_ins, j+2, i-2) ;
  a11 = tex2D(tex_img_ins, j+3, i-2) ;
  //troisieme ligne
  a12 = tex2D(tex_img_ins, j-2, i-1) ;
  a13 = tex2D(tex_img_ins, j-1, i-1) ;
  a14 = tex2D(tex_img_ins, j  , i-1) ;
  a15 = tex2D(tex_img_ins, j+1, i-1) ;
  a16 = tex2D(tex_img_ins, j+2, i-1) ;
  a17 = tex2D(tex_img_ins, j+3, i-1) ;
  //quatrieme ligne
  a18 = tex2D(tex_img_ins, j-2, i  ) ;
  a19 = tex2D(tex_img_ins, j-1, i  ) ;
  a20 = tex2D(tex_img_ins, j  , i  ) ;
  a21 = tex2D(tex_img_ins, j+1, i  ) ;
  a22 = tex2D(tex_img_ins, j+2, i  ) ;
  a23 = tex2D(tex_img_ins, j+3, i  ) ;
  //cinquieme ligne
  a24 = tex2D(tex_img_ins, j-2, i+1) ;
  a25 = tex2D(tex_img_ins, j-1, i+1) ;
  
  //min max aux extremites
  minmax26(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeurs suivante (26)
  a25 = tex2D(tex_img_ins, j , i+1);
  //minmax aux extremites
  minmax25(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (27)
  a25 = tex2D(tex_img_ins, j+1, i+1);
  //minmax aux extremites
  minmax24(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (28)
  a25 = tex2D(tex_img_ins, j+2, i+1);
  //minmax aux extremites
  minmax23(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (29)
  a25 = tex2D(tex_img_ins, j+3, i+1);
  //minmax aux extremites
  minmax22(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (30)
  a25 = tex2D(tex_img_ins, j-2, i+2);
  //minmax aux extremites
  minmax21(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (31)
  a25 = tex2D(tex_img_ins, j-1, i+2);
  //minmax aux extmites
  minmax20(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (32)
  a25 = tex2D(tex_img_ins, j  , i+2);
  //minmax aux extmites
  minmax19(&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (33)
  a25 = tex2D(tex_img_ins, j+1, i+2);
  //minmax aux extmites
  minmax18(&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (34)
  a25 = tex2D(tex_img_ins, j+2, i+2);
  //minmax aux extmites
  minmax17(&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (35)
  a25 = tex2D(tex_img_ins, j+3, i+2);
  //minmax aux extmites
  minmax16(&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (36)
  a25 = tex2D(tex_img_ins, j-2, i+3);
  //minmax aux extmites
  minmax15(&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (37)
  a25 = tex2D(tex_img_ins, j-1, i+3);
  //minmax aux extmites
  
  minmax14(&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  //chargement valeur suivante (38)
  a25 = tex2D(tex_img_ins, j  , i+3);
  //minmax aux extmites
  minmax13(&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (39)
  a25 = tex2D(tex_img_ins, j+1, i+3);
  //minmax aux extmites
  minmax12(&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (40)
  a25 = tex2D(tex_img_ins, j+2, i+3);
  //minmax aux extmites
  minmax11(&a15, &a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (41)
  a25 = tex2D(tex_img_ins, j+3, i+3);
  //minmax aux extmites
  minmax10(&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  // fin des pixels voisins communs deux pixels centraux
  b17=a17; b18=a18; b19=a19; b20=a20; b21=a21; b22=a22; b23=a23; b24=a24; b25=a25;
  
  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_ins, j-3, i+3);
  b25 = tex2D(tex_img_ins, j+4, i+3);
  //minmax aux extremites
  minmax9(&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  minmax9(&b17,&b18,&b19,&b20,&b21,&b22,&b23,&b24,&b25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_ins, j-3, i+2);
  b25 = tex2D(tex_img_ins, j+4, i+2);
  //minmax aux extremites
  minmax8(&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  minmax8(&b18,&b19,&b20,&b21,&b22,&b23,&b24,&b25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_ins, j-3, i+1);
  b25 = tex2D(tex_img_ins, j+4, i+1);
  //minmax aux extremites
  minmax7(&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  minmax7(&b19,&b20,&b21,&b22,&b23,&b24,&b25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_ins, j-3, i);
  b25 = tex2D(tex_img_ins, j+4, i);
  //minmax aux extremites
  minmax6(&a20,&a21,&a22,&a23,&a24,&a25);
  minmax6(&b20,&b21,&b22,&b23,&b24,&b25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_ins, j-3, i-1);
  b25 = tex2D(tex_img_ins, j+4, i-1);
  //minmax aux extremites
  minmax5(&a21,&a22,&a23,&a24,&a25);
  minmax5(&b21,&b22,&b23,&b24,&b25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_ins, j-3, i-2);
  b25 = tex2D(tex_img_ins, j+4, i-2);
  //minmax aux extremites
  minmax4(&a22,&a23,&a24,&a25);
  minmax4(&b22,&b23,&b24,&b25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_ins, j-3, i-3);
  b25 = tex2D(tex_img_ins, j+4, i-3);
  //minmax aux extremites
  minmax3(&a23,&a24,&a25);
  minmax3(&b23,&b24,&b25);
  
  //medians au milieu !
  output[ __mul24(i, j_dim) +j ] = a24 ;
  output[ __mul24(i, j_dim) +j+1 ] = b24 ;
 
}


__global__ void kernel_median7_2pix( unsigned char  *output, int i_dim, int j_dim)
{  
  
  // coordonnees absolues du point
  int j = __mul24(__mul24(blockIdx.x,blockDim.x) + threadIdx.x,2) ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ;
  
  /**************************************************************************
   *             tri(s)
   **************************************************************************/
  int a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,a25 ;
  int b17,b18,b19,b20,b21,b22,b23,b24,b25;															

  /****************************************
   * les 26 premieres valeurs 
   ****************************************/
  //premiere ligne
  a0  = tex2D(tex_img_inc, j-2, i-3) ;
  a1  = tex2D(tex_img_inc, j-1, i-3) ;
  a2  = tex2D(tex_img_inc, j  , i-3) ;
  a3  = tex2D(tex_img_inc, j+1, i-3) ;
  a4  = tex2D(tex_img_inc, j+2, i-3) ;
  a5  = tex2D(tex_img_inc, j+3, i-3) ;
  //deuxieme ligne
  a6  = tex2D(tex_img_inc, j-2, i-2) ;
  a7  = tex2D(tex_img_inc, j-1, i-2) ;
  a8  = tex2D(tex_img_inc, j  , i-2) ;
  a9  = tex2D(tex_img_inc, j+1, i-2) ;
  a10 = tex2D(tex_img_inc, j+2, i-2) ;
  a11 = tex2D(tex_img_inc, j+3, i-2) ;
  //troisieme ligne
  a12 = tex2D(tex_img_inc, j-2, i-1) ;
  a13 = tex2D(tex_img_inc, j-1, i-1) ;
  a14 = tex2D(tex_img_inc, j  , i-1) ;
  a15 = tex2D(tex_img_inc, j+1, i-1) ;
  a16 = tex2D(tex_img_inc, j+2, i-1) ;
  a17 = tex2D(tex_img_inc, j+3, i-1) ;
  //quatrieme ligne
  a18 = tex2D(tex_img_inc, j-2, i  ) ;
  a19 = tex2D(tex_img_inc, j-1, i  ) ;
  a20 = tex2D(tex_img_inc, j  , i  ) ;
  a21 = tex2D(tex_img_inc, j+1, i  ) ;
  a22 = tex2D(tex_img_inc, j+2, i  ) ;
  a23 = tex2D(tex_img_inc, j+3, i  ) ;
  //cinquieme ligne
  a24 = tex2D(tex_img_inc, j-2, i+1) ;
  a25 = tex2D(tex_img_inc, j-1, i+1) ;
  
  //min max aux extremites
  minmax26(&a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeurs suivante (26)
  a25 = tex2D(tex_img_inc, j , i+1);
  //minmax aux extremites
  minmax25(&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (27)
  a25 = tex2D(tex_img_inc, j+1, i+1);
  //minmax aux extremites
  minmax24(&a2,&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (28)
  a25 = tex2D(tex_img_inc, j+2, i+1);
  //minmax aux extremites
  minmax23(&a3,&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (29)
  a25 = tex2D(tex_img_inc, j+3, i+1);
  //minmax aux extremites
  minmax22(&a4,&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (30)
  a25 = tex2D(tex_img_inc, j-2, i+2);
  //minmax aux extremites
  minmax21(&a5,&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (31)
  a25 = tex2D(tex_img_inc, j-1, i+2);
  //minmax aux extmites
  minmax20(&a6,&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (32)
  a25 = tex2D(tex_img_inc, j  , i+2);
  //minmax aux extmites
  minmax19(&a7,&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (33)
  a25 = tex2D(tex_img_inc, j+1, i+2);
  //minmax aux extmites
  minmax18(&a8,&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (34)
  a25 = tex2D(tex_img_inc, j+2, i+2);
  //minmax aux extmites
  minmax17(&a9,&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (35)
  a25 = tex2D(tex_img_inc, j+3, i+2);
  //minmax aux extmites
  minmax16(&a10,&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (36)
  a25 = tex2D(tex_img_inc, j-2, i+3);
  //minmax aux extmites
  minmax15(&a11,&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  //chargement valeur suivante (37)
  a25 = tex2D(tex_img_inc, j-1, i+3);
  //minmax aux extmites
  
  minmax14(&a12,&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  //chargement valeur suivante (38)
  a25 = tex2D(tex_img_inc, j  , i+3);
  //minmax aux extmites
  minmax13(&a13,&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (39)
  a25 = tex2D(tex_img_inc, j+1, i+3);
  //minmax aux extmites
  minmax12(&a14,&a15,&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (40)
  a25 = tex2D(tex_img_inc, j+2, i+3);
  //minmax aux extmites
  minmax11(&a15, &a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);

  //chargement valeur suivante (41)
  a25 = tex2D(tex_img_inc, j+3, i+3);
  //minmax aux extmites
  minmax10(&a16,&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  
  // fin des pixels voisins communs deux pixels centraux
  b17=a17; b18=a18; b19=a19; b20=a20; b21=a21; b22=a22; b23=a23; b24=a24; b25=a25;
  
  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_inc, j-3, i+3);
  b25 = tex2D(tex_img_inc, j+4, i+3);
  //minmax aux extremites
  minmax9(&a17,&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  minmax9(&b17,&b18,&b19,&b20,&b21,&b22,&b23,&b24,&b25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_inc, j-3, i+2);
  b25 = tex2D(tex_img_inc, j+4, i+2);
  //minmax aux extremites
  minmax8(&a18,&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  minmax8(&b18,&b19,&b20,&b21,&b22,&b23,&b24,&b25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_inc, j-3, i+1);
  b25 = tex2D(tex_img_inc, j+4, i+1);
  //minmax aux extremites
  minmax7(&a19,&a20,&a21,&a22,&a23,&a24,&a25);
  minmax7(&b19,&b20,&b21,&b22,&b23,&b24,&b25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_inc, j-3, i);
  b25 = tex2D(tex_img_inc, j+4, i);
  //minmax aux extremites
  minmax6(&a20,&a21,&a22,&a23,&a24,&a25);
  minmax6(&b20,&b21,&b22,&b23,&b24,&b25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_inc, j-3, i-1);
  b25 = tex2D(tex_img_inc, j+4, i-1);
  //minmax aux extremites
  minmax5(&a21,&a22,&a23,&a24,&a25);
  minmax5(&b21,&b22,&b23,&b24,&b25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_inc, j-3, i-2);
  b25 = tex2D(tex_img_inc, j+4, i-2);
  //minmax aux extremites
  minmax4(&a22,&a23,&a24,&a25);
  minmax4(&b22,&b23,&b24,&b25);

  //chargement valeurs suivantes 
  a25 = tex2D(tex_img_inc, j-3, i-3);
  b25 = tex2D(tex_img_inc, j+4, i-3);
  //minmax aux extremites
  minmax3(&a23,&a24,&a25);
  minmax3(&b23,&b24,&b25);
  
  //medians au milieu !
  output[ __mul24(i, j_dim) +j ] = a24 ;
  output[ __mul24(i, j_dim) +j+1 ] = b24 ;
 
}



/*****************************************************************************
 *      median generic shared mem + forgetfull
 *****************************************************************************/
__global__ void kernel_medianForgetRSH( unsigned short * output, int i_dim, int j_dim, int r)
{   
  // coordonnees absolues du point
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ; 
  unsigned short ic, jc ;
  unsigned short Nreg = ((2*r+1)*(2*r+1))/2 + 2 ; 

  int bdimX = blockDim.x; //<<1 ; // pour le cas \E0 2 pixels par thread
  int tidX = threadIdx.x; //<<1 ;
  
  // chargement en smem
  int idrow = threadIdx.y*(blockDim.x+2*r) ;
  
  extern __shared__ int medRoi[];

  // bloc 0 (en haut \E0 gauche) 
  medRoi[ idrow  + threadIdx.x ] = tex2D(tex_img_ins, j-r, i-r) ;
  // bloc 1 (en haut \E0 droite)...
  if ( threadIdx.x < 2*r ) //...ou plutot ce qu'il en manque
	medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-r, i-r ) ;
  // bloc 2 ( en bas \E0 gauche)
  if ( threadIdx.y < 2*r )
	{
	  idrow = (threadIdx.y+blockDim.y)*(blockDim.x+2*r) ;
	  medRoi[ idrow + threadIdx.x ] = tex2D( tex_img_ins, j-r, i+blockDim.y-r ) ;
	  //bloc 4 ( en bas \E0 doite )...
	  if ( threadIdx.x < 2*r ) //...ou ce qu'il en manque
		medRoi[ idrow + threadIdx.x + blockDim.x ] = tex2D( tex_img_ins, j+blockDim.x-r, i+blockDim.y-r ) ;
	}
  __syncthreads() ;

  // remplissage du vecteur de tri minmax
  unsigned short vect[8066] ;
  int Freg=Nreg ;
  for (ic=0; ic<2*r+1; ic++)
	{
	  for (jc=0; jc<2*r+1; jc++)
		{
		  if ( ic*(2*r+1)+jc < Nreg )
			{
			  vect[ ic*(2*r+1)+jc ] = medRoi[ (threadIdx.y+ic)*(bdimX+2*r)+ (tidX+jc) ] ;
			} else
			{
			  minmaxN(vect, Freg--) ;
			  vect[ Nreg-1 ] = medRoi[ (threadIdx.y+ic)*(bdimX+2*r)+ (tidX+jc) ] ;
			}
		}
	}
  minmax3(&vect[Nreg-3], &vect[Nreg-2], &vect[Nreg-1])
  
  //medRoi[ (threadIdx.y+ic)*(bdimX+L-1)+ (tidX+jc) ]
  
	output[ __mul24(i, j_dim) +j ] = vect[ Nreg-2 ];  
}


/*****************************************************************************
 *      median generic shared mem + forgetfull
 *****************************************************************************/
__global__ void kernel_medianForgetR( unsigned short * output, int i_dim, int j_dim, int r)
{   
  // coordonnees absolues du point
  int j = __mul24(blockIdx.x,blockDim.x) + threadIdx.x ; 
  int i = __mul24(blockIdx.y,blockDim.y) + threadIdx.y ; 
  unsigned short ic, jc ;
  unsigned short Nreg = ((2*r+1)*(2*r+1))/2 + 2 ; 

  // remplissage du vecteur de tri minmax
  unsigned short vect[8066] ;
  int Freg=Nreg ;
  for (ic=0; ic<2*r+1; ic++)
	{
	  for (jc=0; jc<2*r+1; jc++)
		{
		  if ( ic*(2*r+1)+jc < Nreg )
			{
			  vect[ ic*(2*r+1)+jc ] = tex2D(tex_img_ins, j-r+jc, i-r+ic) ;
			} else
			{
			  minmaxN(vect, Freg--) ;
			  vect[ Nreg-1 ] = tex2D(tex_img_ins, j-r+jc, i-r+ic) ;
			}
		}
	}
  minmax3(&vect[Nreg-3], &vect[Nreg-2], &vect[Nreg-1])
  
  //medRoi[ (threadIdx.y+ic)*(bdimX+L-1)+ (tidX+jc) ]
  
	output[ __mul24(i, j_dim) +j ] = vect[ Nreg-2 ];  
}



/**********************************************************************************
 *    MEDIAN PSEUDO-SEPARABLE POUR GRANDS NOYAUX
 **********************************************************************************/
__global__ void kernel_medianSepR( unsigned short *output, int i_dim, int j_dim, int r)
{
  
  int idc, val, min, max, inf, egal, sup, mxinf, minsup, estim ;

  //coordonnées ds le bloc
  int ib = threadIdx.y ;
  int jb = threadIdx.x ;
  //int idx_h = __mul24(ib+r,blockDim.x) + jb ;   // index pixel deans shmem (bloc+halo)
  //int offset = __mul24(blockDim.x,r) ;
  
  // coordonnees absolues du point
  int j = __mul24(blockIdx.x,blockDim.x) + jb ; 
  int i = __mul24(blockIdx.y,blockDim.y) + ib ;
  
  extern __shared__ int buff[] ;
  /***********************************************************************************
   *              CHARGEMENT DATA EN SHARED MEM
   ***********************************************************************************/
  for (idc = 0 ; idc < (2*(blockDim.y+r)-1)/blockDim.y; idc++)
	{
	  if (idc*blockDim.y +ib < i_dim)
		buff[ (idc*blockDim.y +ib)*blockDim.x + jb ]    = tex2D(tex_img_ins, j, i-r+ idc*blockDim.y) ;
	}
  
  __syncthreads() ;
  /**********************************************************************************************
   *               TRI VERTICAL par algo TORBEN MOGENSEN
   *          (a little bit slow but saves memory => faster !)
   **********************************************************************************************/
  min = max = buff[ ib*blockDim.x +jb] ;
  
  for (idc= 0 ; idc< 2*r+1 ; idc++ )
	{
	  val = buff[ __mul24(ib+idc, blockDim.x) +jb ] ;
	  if ( val  < min ) min = val ;
	  if ( val  > max ) max = val ;
	}
  
  while (1)
	{  
	  estim = (min+max)/2 ;
	  inf = sup = egal = 0  ;
	  mxinf = min ;
	  minsup= max ;
	  for (idc =0; idc< 2*r+1 ; idc++)
		{
		  val = buff[ __mul24(ib+idc, blockDim.x) +jb ] ;
		  if( val < estim )
			{
			  inf++;
			  if( val > mxinf) mxinf = val ;
			} else if (val > estim)
			{
			  sup++;
			  if( val < minsup) minsup = val ;
			} else egal++ ;
		}
	  if ( (inf <= (r+1))&&(sup <=(r+1)) ) break ;
	  else if (inf>sup) max = mxinf ;
	  else min = minsup ;
	  }
  
  if ( inf >= r+1 ) val = mxinf ;
  else if (inf+egal >= r+1) val = estim ;
  else val = minsup ;
  
  output[ __mul24(j, i_dim)  +i  ] =  val ; 
}


/**
 *
 * correction de staircase
 */
__global__ void kernel_staircase_reduc3(unsigned int * img_out, unsigned int L, unsigned int H)
{
  // coordonnees du point dans le bloc
  //unsigned int iib = threadIdx.x ;
  //unsigned int jib = threadIdx.y ;
  // coordonnees du point dans l'image
  unsigned int y = blockIdx.y*blockDim.y + threadIdx.y;
  unsigned int x = blockIdx.x*blockDim.x + threadIdx.x;

  int a, b, c, d, e, f, g, h, i ;            // gl des voisins 
  float wa, wb, wc, wd, we, wf, wg, wh, wi ; // poids
  float S1, S2, S11, S22, S12, S0, Sx, Sx1, Sx2 ; 
  float c1,c2 ; 

  // chargement des valeurs GL des pixels du voisinage
  a = tex2D(tex_img_estim, x-1, y-1) ;
  b = tex2D(tex_img_estim, x  , y-1) ;
  c = tex2D(tex_img_estim, x+1, y-1) ;
  d = tex2D(tex_img_estim, x-1, y  ) ;
  e = tex2D(tex_img_estim, x  , y  ) ;
  f = tex2D(tex_img_estim, x+1, y  ) ;
  g = tex2D(tex_img_estim, x-1, y+1) ;
  h = tex2D(tex_img_estim, x  , y+1) ;
  i = tex2D(tex_img_estim, x+1, y+1) ;

  wa = tex1D(tex_noyau, abs(a-e)) ;
  wb = tex1D(tex_noyau, abs(b-e)) ;
  wc = tex1D(tex_noyau, abs(c-e)) ;
  wd = tex1D(tex_noyau, abs(d-e)) ;
  we = 1 ;
  wf = tex1D(tex_noyau, abs(f-e)) ;
  wg = tex1D(tex_noyau, abs(g-e)) ;
  wh = tex1D(tex_noyau, abs(h-e)) ;
  wi = tex1D(tex_noyau, abs(i-e)) ;

  
  //construction des elements du systeme lineaire
  S0  = wa+wb+wc+wd+we+wf+wg+wh+wi ; 
  S1  = wc+wf+wi-wa-wd-wg ; 
  S2  = wg+wh+wi-wa-wb-wc ; 
  Sx  = wa*a + wb*b + wc*c + wd*d + we*e + wf*f + wg*g + wh*h + wi*i ;
  Sx1 = wc*c + wf*f + wi*i - wa*a - wd*d - wg*g ;
  Sx2 = wg*g + wh*h + wi*i - wa*a - wb*b - wc*c ;
  S11 = wc+wf+wi+wa+wd+wg ;
  S22 = wg+wh+wi+wa+wb+wc ;
  S12 = wa + wi -wc -wg ;

  c1 = S22*(S11*Sx-S1*Sx1) + S12*(S1*Sx2-S12*Sx) + S2*(S12*Sx1-S11*Sx2) ;
  c2 = S22*(S11*S0-S1*S1) + S12*(S2*S1-S12*S0) + S2*(S1*S12-S2*S11) ;
  img_out[y*L + x] = c1/c2 ;
}