paper_fast_median/code/kernMedianSeparable.cu

   1 __global__ void kernel_medianV_sh( short *output, int i_dim, int j_dim, int r)
   2 {
   3
   4   int idc, val, min, max, inf, egal, sup, mxinf, minsup, estim ;
   5
   6   //coordinates in the block
   7   int ib = threadIdx.y ;
   8   int jb = threadIdx.x ;
   9   int idx_h = __mul24(ib+r,blockDim.x) +jb; // base pixel index
  10   int offset = __mul24(blockDim.x,r) ;
  11
  12   int j = __mul24(blockIdx.x,blockDim.x) + jb ;
  13   int i = __mul24(blockIdx.y,blockDim.y) + ib ;
  14
  15   //      DATA PREFETCHING INTO SHARED MEM
  16   extern __shared__ int buff[] ;
  17   buff[ idx_h ] = tex2D(tex_img_ins, j, i) ;
  18
  19   if (ib < r)
  20         {
  21           buff[ idx_h - offset ]    = tex2D(tex_img_ins, j, i-r) ;
  22         } else
  23         if (ib >= (blockDim.y-r))
  24           {
  25                 buff[ idx_h + offset  ]    = tex2D(tex_img_ins, j, i+r) ;
  26           }
  27
  28   __syncthreads() ;
  29
  30   //      TORBEN MOGENSEN SORTING
  31   min = max = buff[ ib*blockDim.x +jb] ;
  32
  33   for (idc= 0 ; idc< 2*r+1 ; idc++ )
  34         {
  35           val = buff[ __mul24(ib+idc, blockDim.x) +jb ] ;
  36           if ( val  < min ) min = val ;
  37           if ( val  > max ) max = val ;
  38         }
  39
  40   while (1)
  41         {
  42           estim = (min+max)/2 ;
  43           inf = sup = egal = 0  ;
  44           mxinf = min ;
  45           minsup= max ;
  46           for (idc =0; idc< 2*r+1 ; idc++)
  47                 {
  48                   val = buff[ __mul24(ib+idc, blockDim.x) +jb ] ;
  49                   if( val < estim )
  50                         {
  51                           inf++;
  52                           if( val > mxinf) mxinf = val ;
  53                         } else if (val > estim)
  54                         {
  55                           sup++;
  56                           if( val < minsup) minsup = val ;
  57                         } else egal++ ;
  58                 }
  59           if ( (inf <= (r+1))&&(sup <=(r+1)) ) break ;
  60           else if (inf>sup) max = mxinf ;
  61           else min = minsup ;
  62           }
  63
  64   if ( inf >= r+1 ) val = mxinf ;
  65   else if (inf+egal >= r+1) val = estim ;
  66   else val = minsup ;
  67
  68   output[ __mul24(j, i_dim)  +i  ] =  val ;
  69 }