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

/*********************************************************************
 *    MACROS pour tri min-max utilisant la fonction de tri s(a,b) 
 *********************************************************************/
#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);
#define minmax4(a, b, c, d) s(a, b); s(c, d); s(a, c); s(b, d);
#define minmax5(a, b, c, d, e) s(a, b); s(c, d); min3(a, c, e); max3(b, d, e);
#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);

__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 ;
  
  // registres
  int a0, a1, a2, a3, a4, a5 ;
  int b0, b1, b2, b3, b4, b5 ;

  // les six premieres valeurs 
  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) ;
  
  //identification des extrema
  minmax6(&a0, &a1, &a2, &a3, &a4, &a5) ;

  // séparation des processus   
  b0=a0; b1=a1; b2=a2; b3=a3; b4=a4; b5=a5;
  
  // les valeurs suivantes en fin de liste 
  a5 = tex2D(tex_img_ins, j-1, i  ) ;
  b5 = tex2D(tex_img_ins, j+2, i  ) ;
  // identification des extrema
  minmax5(&a1, &a2, &a3, &a4, &a5) ;
  minmax5(&b1, &b2, &b3, &b4, &b5) ;

  // les 2 dernières étapes de la sélection
  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) ;
  
  // les medians sont en position milieu
  // on sort les 2 valeurs à la suite 
  output[ __mul24(i, j_dim) +j ] = a4 ;
  output[ __mul24(i, j_dim) +j+1 ] = b4 ;
 
}