def get_kernel_source_template():
    """
    Returns the OpenCL 2D median kernel source based on `Viola et al. 2003
    <http://dk.migrid.org/vgrid/eScience/Tools/Python/litterature/
    cphhpctoolbox/cphhpc/image/filter/Viola-2003.pdf>`_ section 4.2.1.
    """
    return """

#ifdef DATA_TYPE_UINT16
    __kernel void median2d(__global unsigned short *result, 
                            __global unsigned short *input) {

                            
       __local unsigned short shared_window[SHARED_WINDOW_HEIGHT][SHARED_WINDOW_WIDTH];
       
#elif defined DATA_TYPE_UINT32
    __kernel void median2d(__global unsigned int *result, 
                            __global unsigned int *input) {

                            
       __local unsigned int shared_window[SHARED_WINDOW_HEIGHT][SHARED_WINDOW_WIDTH];
#endif


       int x = get_global_id(0);
       int y = get_global_id(1);
       
       int tx = get_local_id(0);
       int ty = get_local_id(1);
        
       int tx_top = tx-FILTER_HALF_X_SIZE;
       int tx_bottom = tx+FILTER_HALF_X_SIZE;
        
       int ty_top = ty-FILTER_HALF_Y_SIZE;
       int ty_bottom = ty+FILTER_HALF_Y_SIZE;
       
       /* 
       * zero pad block borders 
       */
       SHARED_WINDOW(ty, tx_top) = 0;
       SHARED_WINDOW(ty, tx_bottom) = 0;
       SHARED_WINDOW(ty_top, tx) = 0;
       SHARED_WINDOW(ty_top, tx_top) = 0;
       SHARED_WINDOW(ty_top, tx_bottom) = 0;
       SHARED_WINDOW(ty_bottom, tx) = 0;
       SHARED_WINDOW(ty_bottom, tx_top) = 0;
       SHARED_WINDOW(ty_bottom, tx_bottom) = 0;
      
       /*
       * As above zero pad pads non-border indexes (center threads)
       * we need to sync before proceed fetching data
       * NOTE: this is more efficient than checking if we are at a border
       */
       barrier(CLK_LOCAL_MEM_FENCE);
 
       /*
       * guards: is at boundary and still more image?
       */
       bool is_x_top = (tx_top<0);
       bool is_y_top = (ty_top<0);
       bool is_x_bot = (tx_bottom>WORK_GROUP_W-1);
       bool is_y_bot = (ty_bottom>WORK_GROUP_H-1);
       
       int x_top = x-FILTER_HALF_X_SIZE;
       int x_bottom = x+FILTER_HALF_X_SIZE;
       
       int y_top = y-FILTER_HALF_Y_SIZE;
       int y_bottom = y+FILTER_HALF_Y_SIZE;
       
       is_x_top = (x_top >= 0);
       is_y_top = (y_top >= 0);
       is_x_bot = (x_bottom<DATA_W);
       is_y_bot = (y_bottom<DATA_H);
       
       /*
       * Each thread pulls from image
       */
      
       // Self 
       SHARED_WINDOW(ty,tx) = INPUT(y,x);
       
       // Borders
       if (is_x_top) {
          SHARED_WINDOW(ty, tx_top) = INPUT(y, x_top);
       }
       
       if (is_x_bot) {
          SHARED_WINDOW(ty, tx_bottom) = INPUT(y, x_bottom);
       }

       if (is_y_top) {
          SHARED_WINDOW(ty_top, tx) = INPUT(y_top, x);
          if (is_x_top) {
             SHARED_WINDOW(ty_top, tx_top) = INPUT(y_top, x_top);
          }
          if (is_x_bot) {
             SHARED_WINDOW(ty_top, tx_bottom) = INPUT(y_top, x_bottom);
          }
       }
       if (is_y_bot) {
          SHARED_WINDOW(ty_bottom, tx) = INPUT(y_bottom, x);
          
          if (is_x_top) {
             SHARED_WINDOW(ty_bottom, tx_top) = INPUT(y_bottom, x_top);
          }
         
          if (is_x_bot) {
             SHARED_WINDOW(ty_bottom, tx_bottom) = INPUT(y_bottom, x_bottom);
          }
       }
      
       /*
        * Wait for threads to finish data fetch
       */ 
       barrier(CLK_LOCAL_MEM_FENCE);

       /*
        * With UINT32 float is _NOT_precise enough
       */
#ifdef DATA_TYPE_UINT16
       float minval = 0.0f;
       float maxval = convert_float(INIT_MAXVAL);
       float pivot = maxval/2.0f;
       
#elif defined DATA_TYPE_UINT32
       double minval = 0.0;
       double maxval = convert_double(INIT_MAXVAL);
       double pivot = maxval/2.0;
#endif

       unsigned short count = 0;
       
       short pos_x;
       short start_x = tx-FILTER_HALF_X_SIZE;
       short end_x = tx+FILTER_HALF_X_SIZE_PLUS_ONE;
      
       short pos_y;
       short start_y = ty-FILTER_HALF_Y_SIZE;
       short end_y = ty+FILTER_HALF_Y_SIZE_PLUS_ONE;
       
       for (unsigned short i=0; i<ITERATIONS; i++) {
          count = 0;
       
          for (pos_y=start_y; pos_y<end_y; pos_y++) {
             for (pos_x=start_x; pos_x<end_x; pos_x++) {
                count = (SHARED_WINDOW(pos_y, pos_x) > pivot)?count+1:count;
             }
          }

          maxval = (count < MEDIAN_PIVOT_SIZE)?floor(pivot):maxval;
          minval = (count >= MEDIAN_PIVOT_SIZE)?ceil(pivot):minval;
          pivot = (minval + maxval)/2.0f;

       }
       RESULT(y,x) = ceil(pivot);

    }
    """
