modules/ocl/src/opencl/moments.cl

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                           License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2010-2012, Multicoreware, Inc., all rights reserved.
  14 // Copyright (C) 2010-2012, Advanced Micro Devices, Inc., all rights reserved.
  15 // Third party copyrights are property of their respective owners.
  16 //
  17 // @Authors
  18 //    Sen Liu, swjtuls1987@126.com
  19 //
  20 // Redistribution and use in source and binary forms, with or without modification,
  21 // are permitted provided that the following conditions are met:
  22 //
  23 //   * Redistribution's of source code must retain the above copyright notice,
  24 //     this list of conditions and the following disclaimer.
  25 //
  26 //   * Redistribution's in binary form must reproduce the above copyright notice,
  27 //     this list of conditions and the following disclaimer in the documentation
  28 //     and/or other oclMaterials provided with the distribution.
  29 //
  30 //   * The name of the copyright holders may not be used to endorse or promote products
  31 //     derived from this software without specific prior written permission.
  32 //
  33 // This software is provided by the copyright holders and contributors as is and
  34 // any express or implied warranties, including, but not limited to, the implied
  35 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  36 // In no event shall the Intel Corporation or contributors be liable for any direct,
  37 // indirect, incidental, special, exemplary, or consequential damages
  38 // (including, but not limited to, procurement of substitute goods or services;
  39 // loss of use, data, or profits; or business interruption) however caused
  40 // and on any theory of liability, whether in contract, strict liability,
  41 // or tort (including negligence or otherwise) arising in any way out of
  42 // the use of this software, even if advised of the possibility of such damage.
  43 //
  44 //M*/
  45
  46 #if defined (DOUBLE_SUPPORT)
  47
  48 #ifdef cl_khr_fp64
  49 #pragma OPENCL EXTENSION cl_khr_fp64:enable
  50 #elif defined (cl_amd_fp64)
  51 #pragma OPENCL EXTENSION cl_amd_fp64:enable
  52 #endif
  53 typedef double T;
  54 typedef double F;
  55 typedef double4 F4;
  56 #define convert_F4 convert_double4
  57
  58 #else
  59 typedef float F;
  60 typedef float4 F4;
  61 typedef long T;
  62 #define convert_F4 convert_float4
  63 #endif
  64
  65 #define DST_ROW_00     0
  66 #define DST_ROW_10     1
  67 #define DST_ROW_01     2
  68 #define DST_ROW_20     3
  69 #define DST_ROW_11     4
  70 #define DST_ROW_02     5
  71 #define DST_ROW_30     6
  72 #define DST_ROW_21     7
  73 #define DST_ROW_12     8
  74 #define DST_ROW_03     9
  75
  76 __kernel void icvContourMoments(int contour_total,
  77                                 __global float* reader_oclmat_data,
  78                                 __global T* dst_a,
  79                                 int dst_step)
  80 {
  81     T xi_1, yi_1, xi_12, yi_12, xi, yi, xi2, yi2, dxy, xii_1, yii_1;
  82     int idx = get_global_id(0);
  83
  84     if (idx < 0 || idx >= contour_total)
  85         return;
  86
  87     xi_1 = (T)(*(reader_oclmat_data + (get_global_id(0) << 1)));
  88     yi_1 = (T)(*(reader_oclmat_data + (get_global_id(0) << 1) + 1));
  89     xi_12 = xi_1 * xi_1;
  90     yi_12 = yi_1 * yi_1;
  91
  92     if(idx == contour_total - 1)
  93     {
  94         xi = (T)(*(reader_oclmat_data));
  95         yi = (T)(*(reader_oclmat_data + 1));
  96     }
  97     else
  98     {
  99         xi = (T)(*(reader_oclmat_data + (idx + 1) * 2));
 100         yi = (T)(*(reader_oclmat_data + (idx + 1) * 2 + 1));
 101     }
 102
 103     xi2 = xi * xi;
 104     yi2 = yi * yi;
 105     dxy = xi_1 * yi - xi * yi_1;
 106     xii_1 = xi_1 + xi;
 107     yii_1 = yi_1 + yi;
 108
 109     dst_step /= sizeof(T);
 110     *( dst_a + DST_ROW_00 * dst_step + idx) = dxy;
 111     *( dst_a + DST_ROW_10 * dst_step + idx) = dxy * xii_1;
 112     *( dst_a + DST_ROW_01 * dst_step + idx) = dxy * yii_1;
 113     *( dst_a + DST_ROW_20 * dst_step + idx) = dxy * (xi_1 * xii_1 + xi2);
 114     *( dst_a + DST_ROW_11 * dst_step + idx) = dxy * (xi_1 * (yii_1 + yi_1) + xi * (yii_1 + yi));
 115     *( dst_a + DST_ROW_02 * dst_step + idx) = dxy * (yi_1 * yii_1 + yi2);
 116     *( dst_a + DST_ROW_30 * dst_step + idx) = dxy * xii_1 * (xi_12 + xi2);
 117     *( dst_a + DST_ROW_03 * dst_step + idx) = dxy * yii_1 * (yi_12 + yi2);
 118     *( dst_a + DST_ROW_21 * dst_step + idx) =
 119         dxy * (xi_12 * (3 * yi_1 + yi) + 2 * xi * xi_1 * yii_1 +
 120                xi2 * (yi_1 + 3 * yi));
 121     *( dst_a + DST_ROW_12 * dst_step + idx) =
 122         dxy * (yi_12 * (3 * xi_1 + xi) + 2 * yi * yi_1 * xii_1 +
 123                yi2 * (xi_1 + 3 * xi));
 124 }
 125
 126 __kernel void dst_sum(int src_rows, int src_cols, int tile_height, int tile_width, int TILE_SIZE,
 127                       __global F* sum, __global F* dst_m, int dst_step)
 128 {
 129     int gidy = get_global_id(0);
 130     int gidx = get_global_id(1);
 131     int block_y = src_rows/tile_height;
 132     int block_x = src_cols/tile_width;
 133     int block_num;
 134
 135     if(src_rows > TILE_SIZE && src_rows % TILE_SIZE != 0)
 136         block_y ++;
 137     if(src_cols > TILE_SIZE && src_cols % TILE_SIZE != 0)
 138         block_x ++;
 139     block_num = block_y * block_x;
 140     __local F dst_sum[10][128];
 141     if(gidy<128-block_num)
 142         for(int i=0; i<10; i++)
 143             dst_sum[i][gidy+block_num]=0;
 144     barrier(CLK_LOCAL_MEM_FENCE);
 145
 146     dst_step /= sizeof(F);
 147     if(gidy<block_num)
 148     {
 149         dst_sum[0][gidy] = *(dst_m + mad24(DST_ROW_00 * block_y, dst_step, gidy));
 150         dst_sum[1][gidy] = *(dst_m + mad24(DST_ROW_10 * block_y, dst_step, gidy));
 151         dst_sum[2][gidy] = *(dst_m + mad24(DST_ROW_01 * block_y, dst_step, gidy));
 152         dst_sum[3][gidy] = *(dst_m + mad24(DST_ROW_20 * block_y, dst_step, gidy));
 153         dst_sum[4][gidy] = *(dst_m + mad24(DST_ROW_11 * block_y, dst_step, gidy));
 154         dst_sum[5][gidy] = *(dst_m + mad24(DST_ROW_02 * block_y, dst_step, gidy));
 155         dst_sum[6][gidy] = *(dst_m + mad24(DST_ROW_30 * block_y, dst_step, gidy));
 156         dst_sum[7][gidy] = *(dst_m + mad24(DST_ROW_21 * block_y, dst_step, gidy));
 157         dst_sum[8][gidy] = *(dst_m + mad24(DST_ROW_12 * block_y, dst_step, gidy));
 158         dst_sum[9][gidy] = *(dst_m + mad24(DST_ROW_03 * block_y, dst_step, gidy));
 159     }
 160     barrier(CLK_LOCAL_MEM_FENCE);
 161     for(int lsize=64; lsize>0; lsize>>=1)
 162     {
 163         if(gidy<lsize)
 164         {
 165             int lsize2 = gidy + lsize;
 166             for(int i=0; i<10; i++)
 167                 dst_sum[i][gidy] += dst_sum[i][lsize2];
 168         }
 169         barrier(CLK_LOCAL_MEM_FENCE);
 170     }
 171     if(gidy==0)
 172         for(int i=0; i<10; i++)
 173             sum[i] = dst_sum[i][0];
 174 }
 175
 176 __kernel void CvMoments_D0(__global uchar16* src_data, int src_rows, int src_cols, int src_step, int tileSize_width, int tileSize_height,
 177                            __global F* dst_m,
 178                            int dst_cols, int dst_step, int blocky,
 179                            int type, int depth, int cn, int coi, int binary, int TILE_SIZE)
 180 {
 181     uchar tmp_coi[16]; // get the coi data
 182     uchar16 tmp[16];
 183     int VLEN_C = 16;  // vector length of uchar
 184
 185     int gidy = get_global_id(0);
 186     int gidx = get_global_id(1);
 187     int wgidy = get_group_id(0);
 188     int wgidx = get_group_id(1);
 189     int lidy = get_local_id(0);
 190     int lidx = get_local_id(1);
 191     int y = wgidy*TILE_SIZE; // vector length of uchar
 192     int x = wgidx*TILE_SIZE;  // vector length of uchar
 193     int kcn = (cn==2)?2:4;
 194     int rstep = min(src_step, TILE_SIZE);
 195     tileSize_height = min(TILE_SIZE, src_rows - y);
 196     tileSize_width = min(TILE_SIZE, src_cols - x);
 197
 198     if( tileSize_width < TILE_SIZE )
 199         for(int i = tileSize_width; i < rstep; i++ )
 200             *((__global uchar*)src_data+(y+lidy)*src_step+x+i) = 0;
 201     if( coi > 0 )       //channel of interest
 202         for(int i = 0; i < tileSize_width; i += VLEN_C)
 203         {
 204             for(int j=0; j<VLEN_C; j++)
 205                 tmp_coi[j] = *((__global uchar*)src_data+(y+lidy)*src_step+(x+i+j)*kcn+coi-1);
 206             tmp[i/VLEN_C] = (uchar16)(tmp_coi[0],tmp_coi[1],tmp_coi[2],tmp_coi[3],tmp_coi[4],tmp_coi[5],tmp_coi[6],tmp_coi[7],
 207                                       tmp_coi[8],tmp_coi[9],tmp_coi[10],tmp_coi[11],tmp_coi[12],tmp_coi[13],tmp_coi[14],tmp_coi[15]);
 208         }
 209     else
 210         for(int i=0; i < tileSize_width; i+=VLEN_C)
 211             tmp[i/VLEN_C] = *(src_data+(y+lidy)*src_step/VLEN_C+(x+i)/VLEN_C);
 212     uchar16 zero = (uchar16)(0);
 213     uchar16 full = (uchar16)(255);
 214     if( binary )
 215         for(int i=0; i < tileSize_width; i+=VLEN_C)
 216             tmp[i/VLEN_C] = (tmp[i/VLEN_C]!=zero)?full:zero;
 217     F mom[10];
 218     __local int m[10][128];
 219     if(lidy == 0)
 220         for(int i=0; i<10; i++)
 221             for(int j=0; j<128; j++)
 222                 m[i][j]=0;
 223     barrier(CLK_LOCAL_MEM_FENCE);
 224     int lm[10] = {0};
 225     int16 x0 = (int16)(0);
 226     int16 x1 = (int16)(0);
 227     int16 x2 = (int16)(0);
 228     int16 x3 = (int16)(0);
 229     for( int xt = 0 ; xt < tileSize_width; xt+=(VLEN_C) )
 230     {
 231         int16 v_xt = (int16)(xt, xt+1, xt+2, xt+3, xt+4, xt+5, xt+6, xt+7, xt+8, xt+9, xt+10, xt+11, xt+12, xt+13, xt+14, xt+15);
 232         int16 p = convert_int16(tmp[xt/VLEN_C]);
 233         int16 xp = v_xt * p, xxp = xp *v_xt;
 234         x0 += p;
 235         x1 += xp;
 236         x2 += xxp;
 237         x3 += xxp * v_xt;
 238     }
 239     x0.s0 += x0.s1 + x0.s2 + x0.s3 + x0.s4 + x0.s5 + x0.s6 + x0.s7 + x0.s8 + x0.s9 + x0.sa + x0.sb + x0.sc + x0.sd + x0.se + x0.sf;
 240     x1.s0 += x1.s1 + x1.s2 + x1.s3 + x1.s4 + x1.s5 + x1.s6 + x1.s7 + x1.s8 + x1.s9 + x1.sa + x1.sb + x1.sc + x1.sd + x1.se + x1.sf;
 241     x2.s0 += x2.s1 + x2.s2 + x2.s3 + x2.s4 + x2.s5 + x2.s6 + x2.s7 + x2.s8 + x2.s9 + x2.sa + x2.sb + x2.sc + x2.sd + x2.se + x2.sf;
 242     x3.s0 += x3.s1 + x3.s2 + x3.s3 + x3.s4 + x3.s5 + x3.s6 + x3.s7 + x3.s8 + x3.s9 + x3.sa + x3.sb + x3.sc + x3.sd + x3.se + x3.sf;
 243     int py = lidy * ((int)x0.s0);
 244     int sy = lidy*lidy;
 245     int bheight = min(tileSize_height, TILE_SIZE/2);
 246     if(bheight >= TILE_SIZE/2&&lidy > bheight-1&&lidy < tileSize_height)
 247     {
 248         m[9][lidy-bheight] = ((int)py) * sy;  // m03
 249         m[8][lidy-bheight] = ((int)x1.s0) * sy;  // m12
 250         m[7][lidy-bheight] = ((int)x2.s0) * lidy;  // m21
 251         m[6][lidy-bheight] = x3.s0;             // m30
 252         m[5][lidy-bheight] = x0.s0 * sy;        // m02
 253         m[4][lidy-bheight] = x1.s0 * lidy;         // m11
 254         m[3][lidy-bheight] = x2.s0;             // m20
 255         m[2][lidy-bheight] = py;             // m01
 256         m[1][lidy-bheight] = x1.s0;             // m10
 257         m[0][lidy-bheight] = x0.s0;             // m00
 258     }
 259     else if(lidy < bheight)
 260     {
 261         lm[9] = ((int)py) * sy;  // m03
 262         lm[8] = ((int)x1.s0) * sy;  // m12
 263         lm[7] = ((int)x2.s0) * lidy;  // m21
 264         lm[6] = x3.s0;             // m30
 265         lm[5] = x0.s0 * sy;        // m02
 266         lm[4] = x1.s0 * lidy;         // m11
 267         lm[3] = x2.s0;             // m20
 268         lm[2] = py;             // m01
 269         lm[1] = x1.s0;             // m10
 270         lm[0] = x0.s0;             // m00
 271     }
 272     barrier(CLK_LOCAL_MEM_FENCE);
 273     for( int j = bheight; j >= 1; j = j/2 )
 274     {
 275         if(lidy < j)
 276             for( int i = 0; i < 10; i++ )
 277                 lm[i] = lm[i] + m[i][lidy];
 278         barrier(CLK_LOCAL_MEM_FENCE);
 279         if(lidy >= j/2&&lidy < j)
 280             for( int i = 0; i < 10; i++ )
 281                 m[i][lidy-j/2] = lm[i];
 282         barrier(CLK_LOCAL_MEM_FENCE);
 283     }
 284     if(lidy == 0&&lidx == 0)
 285     {
 286         for( int mt = 0; mt < 10; mt++ )
 287             mom[mt] = (F)lm[mt];
 288         if(binary)
 289         {
 290             F s = 1./255;
 291             for( int mt = 0; mt < 10; mt++ )
 292                 mom[mt] *= s;
 293         }
 294         F xm = x * mom[0], ym = y * mom[0];
 295
 296         // accumulate moments computed in each tile
 297         dst_step /= sizeof(F);
 298
 299         // + m00 ( = m00' )
 300         *(dst_m + mad24(DST_ROW_00 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[0];
 301
 302         // + m10 ( = m10' + x*m00' )
 303         *(dst_m + mad24(DST_ROW_10 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[1] + xm;
 304
 305         // + m01 ( = m01' + y*m00' )
 306         *(dst_m + mad24(DST_ROW_01 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[2] + ym;
 307
 308         // + m20 ( = m20' + 2*x*m10' + x*x*m00' )
 309         *(dst_m + mad24(DST_ROW_20 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[3] + x * (mom[1] * 2 + xm);
 310
 311         // + m11 ( = m11' + x*m01' + y*m10' + x*y*m00' )
 312         *(dst_m + mad24(DST_ROW_11 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[4] + x * (mom[2] + ym) + y * mom[1];
 313
 314         // + m02 ( = m02' + 2*y*m01' + y*y*m00' )
 315         *(dst_m + mad24(DST_ROW_02 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[5] + y * (mom[2] * 2 + ym);
 316
 317         // + m30 ( = m30' + 3*x*m20' + 3*x*x*m10' + x*x*x*m00' )
 318         *(dst_m + mad24(DST_ROW_30 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[6] + x * (3. * mom[3] + x * (3. * mom[1] + xm));
 319
 320         // + m21 ( = m21' + x*(2*m11' + 2*y*m10' + x*m01' + x*y*m00') + y*m20')
 321         *(dst_m + mad24(DST_ROW_21 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[7] + x * (2 * (mom[4] + y * mom[1]) + x * (mom[2] + ym)) + y * mom[3];
 322
 323         // + m12 ( = m12' + y*(2*m11' + 2*x*m01' + y*m10' + x*y*m00') + x*m02')
 324         *(dst_m + mad24(DST_ROW_12 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[8] + y * (2 * (mom[4] + x * mom[2]) + y * (mom[1] + xm)) + x * mom[5];
 325
 326         // + m03 ( = m03' + 3*y*m02' + 3*y*y*m01' + y*y*y*m00' )
 327         *(dst_m + mad24(DST_ROW_03 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[9] + y * (3. * mom[5] + y * (3. * mom[2] + ym));
 328     }
 329 }
 330
 331 __kernel void CvMoments_D2(__global ushort8* src_data, int src_rows, int src_cols, int src_step, int tileSize_width, int tileSize_height,
 332                            __global F* dst_m,
 333                            int dst_cols, int dst_step, int blocky,
 334                            int type, int depth, int cn, int coi, int binary, const int TILE_SIZE)
 335 {
 336     ushort tmp_coi[8]; // get the coi data
 337     ushort8 tmp[32];
 338     int VLEN_US = 8; // vector length of ushort
 339     int gidy = get_global_id(0);
 340     int gidx = get_global_id(1);
 341     int wgidy = get_group_id(0);
 342     int wgidx = get_group_id(1);
 343     int lidy = get_local_id(0);
 344     int lidx = get_local_id(1);
 345     int y = wgidy*TILE_SIZE;  // real Y index of pixel
 346     int x = wgidx*TILE_SIZE;  // real X index of pixel
 347     int kcn = (cn==2)?2:4;
 348     int rstep = min(src_step/2, TILE_SIZE);
 349     tileSize_height = min(TILE_SIZE, src_rows - y);
 350     tileSize_width = min(TILE_SIZE, src_cols -x);
 351     if(src_cols > TILE_SIZE && tileSize_width < TILE_SIZE)
 352         for(int i=tileSize_width; i < rstep; i++ )
 353             *((__global ushort*)src_data+(y+lidy)*src_step/2+x+i) = 0;
 354     if( coi > 0 )
 355         for(int i=0; i < tileSize_width; i+=VLEN_US)
 356         {
 357             for(int j=0; j<VLEN_US; j++)
 358                 tmp_coi[j] = *((__global ushort*)src_data+(y+lidy)*(int)src_step/2+(x+i+j)*kcn+coi-1);
 359             tmp[i/VLEN_US] = (ushort8)(tmp_coi[0],tmp_coi[1],tmp_coi[2],tmp_coi[3],tmp_coi[4],tmp_coi[5],tmp_coi[6],tmp_coi[7]);
 360         }
 361     else
 362         for(int i=0; i < tileSize_width; i+=VLEN_US)
 363             tmp[i/VLEN_US] = *(src_data+(y+lidy)*src_step/(2*VLEN_US)+(x+i)/VLEN_US);
 364     ushort8 zero = (ushort8)(0);
 365     ushort8 full = (ushort8)(255);
 366     if( binary )
 367         for(int i=0; i < tileSize_width; i+=VLEN_US)
 368             tmp[i/VLEN_US] = (tmp[i/VLEN_US]!=zero)?full:zero;
 369     F mom[10];
 370     __local long m[10][128];
 371     if(lidy == 0)
 372         for(int i=0; i<10; i++)
 373             for(int j=0; j<128; j++)
 374                 m[i][j]=0;
 375     barrier(CLK_LOCAL_MEM_FENCE);
 376     long lm[10] = {0};
 377     int8 x0 = (int8)(0);
 378     int8 x1 = (int8)(0);
 379     int8 x2 = (int8)(0);
 380     long8 x3 = (long8)(0);
 381     for( int xt = 0 ; xt < tileSize_width; xt+=(VLEN_US) )
 382     {
 383         int8 v_xt = (int8)(xt, xt+1, xt+2, xt+3, xt+4, xt+5, xt+6, xt+7);
 384         int8 p = convert_int8(tmp[xt/VLEN_US]);
 385         int8 xp = v_xt * p, xxp = xp * v_xt;
 386         x0 += p;
 387         x1 += xp;
 388         x2 += xxp;
 389         x3 += convert_long8(xxp) *convert_long8(v_xt);
 390     }
 391     x0.s0 += x0.s1 + x0.s2 + x0.s3 + x0.s4 + x0.s5 + x0.s6 + x0.s7;
 392     x1.s0 += x1.s1 + x1.s2 + x1.s3 + x1.s4 + x1.s5 + x1.s6 + x1.s7;
 393     x2.s0 += x2.s1 + x2.s2 + x2.s3 + x2.s4 + x2.s5 + x2.s6 + x2.s7;
 394     x3.s0 += x3.s1 + x3.s2 + x3.s3 + x3.s4 + x3.s5 + x3.s6 + x3.s7;
 395
 396     int py = lidy * x0.s0, sy = lidy*lidy;
 397     int bheight = min(tileSize_height, TILE_SIZE/2);
 398     if(bheight >= TILE_SIZE/2&&lidy > bheight-1&&lidy < tileSize_height)
 399     {
 400         m[9][lidy-bheight] = ((long)py) * sy;  // m03
 401         m[8][lidy-bheight] = ((long)x1.s0) * sy;  // m12
 402         m[7][lidy-bheight] = ((long)x2.s0) * lidy;  // m21
 403         m[6][lidy-bheight] = x3.s0;             // m30
 404         m[5][lidy-bheight] = x0.s0 * sy;        // m02
 405         m[4][lidy-bheight] = x1.s0 * lidy;         // m11
 406         m[3][lidy-bheight] = x2.s0;             // m20
 407         m[2][lidy-bheight] = py;             // m01
 408         m[1][lidy-bheight] = x1.s0;             // m10
 409         m[0][lidy-bheight] = x0.s0;             // m00
 410     }
 411     else if(lidy < bheight)
 412     {
 413         lm[9] = ((long)py) * sy;  // m03
 414         lm[8] = ((long)x1.s0) * sy;  // m12
 415         lm[7] = ((long)x2.s0) * lidy;  // m21
 416         lm[6] = x3.s0;             // m30
 417         lm[5] = x0.s0 * sy;        // m02
 418         lm[4] = x1.s0 * lidy;         // m11
 419         lm[3] = x2.s0;             // m20
 420         lm[2] = py;             // m01
 421         lm[1] = x1.s0;             // m10
 422         lm[0] = x0.s0;             // m00
 423     }
 424     barrier(CLK_LOCAL_MEM_FENCE);
 425     for( int j = TILE_SIZE/2; j >= 1; j = j/2 )
 426     {
 427         if(lidy < j)
 428             for( int i = 0; i < 10; i++ )
 429                 lm[i] = lm[i] + m[i][lidy];
 430         barrier(CLK_LOCAL_MEM_FENCE);
 431         if(lidy >= j/2&&lidy < j)
 432             for( int i = 0; i < 10; i++ )
 433                 m[i][lidy-j/2] = lm[i];
 434         barrier(CLK_LOCAL_MEM_FENCE);
 435     }
 436     if(lidy == 0&&lidx == 0)
 437     {
 438         for(int mt = 0; mt < 10; mt++ )
 439             mom[mt] = (F)lm[mt];
 440
 441         if(binary)
 442         {
 443             F s = 1./255;
 444             for( int mt = 0; mt < 10; mt++ )
 445                 mom[mt] *= s;
 446         }
 447
 448         F xm = x  *mom[0], ym = y * mom[0];
 449
 450         // accumulate moments computed in each tile
 451         dst_step /= sizeof(F);
 452
 453         // + m00 ( = m00' )
 454         *(dst_m + mad24(DST_ROW_00 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[0];
 455
 456         // + m10 ( = m10' + x*m00' )
 457         *(dst_m + mad24(DST_ROW_10 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[1] + xm;
 458
 459         // + m01 ( = m01' + y*m00' )
 460         *(dst_m + mad24(DST_ROW_01 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[2] + ym;
 461
 462         // + m20 ( = m20' + 2*x*m10' + x*x*m00' )
 463         *(dst_m + mad24(DST_ROW_20 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[3] + x * (mom[1] * 2 + xm);
 464
 465         // + m11 ( = m11' + x*m01' + y*m10' + x*y*m00' )
 466         *(dst_m + mad24(DST_ROW_11 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[4] + x * (mom[2] + ym) + y * mom[1];
 467
 468         // + m02 ( = m02' + 2*y*m01' + y*y*m00' )
 469         *(dst_m + mad24(DST_ROW_02 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[5] + y * (mom[2] * 2 + ym);
 470
 471         // + m30 ( = m30' + 3*x*m20' + 3*x*x*m10' + x*x*x*m00' )
 472         *(dst_m + mad24(DST_ROW_30 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[6] + x * (3. * mom[3] + x * (3. * mom[1] + xm));
 473
 474         // + m21 ( = m21' + x*(2*m11' + 2*y*m10' + x*m01' + x*y*m00') + y*m20')
 475         *(dst_m + mad24(DST_ROW_21 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[7] + x * (2 * (mom[4] + y * mom[1]) + x * (mom[2] + ym)) + y * mom[3];
 476
 477         // + m12 ( = m12' + y*(2*m11' + 2*x*m01' + y*m10' + x*y*m00') + x*m02')
 478         *(dst_m + mad24(DST_ROW_12 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[8] + y * (2 * (mom[4] + x * mom[2]) + y * (mom[1] + xm)) + x * mom[5];
 479
 480         // + m03 ( = m03' + 3*y*m02' + 3*y*y*m01' + y*y*y*m00' )
 481         *(dst_m + mad24(DST_ROW_03 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[9] + y * (3. * mom[5] + y * (3. * mom[2] + ym));
 482     }
 483 }
 484
 485 __kernel void CvMoments_D3(__global short8* src_data, int src_rows, int src_cols, int src_step, int tileSize_width, int tileSize_height,
 486                            __global F* dst_m,
 487                            int dst_cols, int dst_step, int blocky,
 488                            int type, int depth, int cn, int coi, int binary, const int TILE_SIZE)
 489 {
 490     short tmp_coi[8]; // get the coi data
 491     short8 tmp[32];
 492     int VLEN_S =8; // vector length of short
 493     int gidy = get_global_id(0);
 494     int gidx = get_global_id(1);
 495     int wgidy = get_group_id(0);
 496     int wgidx = get_group_id(1);
 497     int lidy = get_local_id(0);
 498     int lidx = get_local_id(1);
 499     int y = wgidy*TILE_SIZE;  // real Y index of pixel
 500     int x = wgidx*TILE_SIZE;  // real X index of pixel
 501     int kcn = (cn==2)?2:4;
 502     int rstep = min(src_step/2, TILE_SIZE);
 503     tileSize_height = min(TILE_SIZE, src_rows - y);
 504     tileSize_width = min(TILE_SIZE, src_cols -x);
 505     if(tileSize_width < TILE_SIZE)
 506         for(int i = tileSize_width; i < rstep; i++ )
 507             *((__global short*)src_data+(y+lidy)*src_step/2+x+i) = 0;
 508     if( coi > 0 )
 509         for(int i=0; i < tileSize_width; i+=VLEN_S)
 510         {
 511             for(int j=0; j<VLEN_S; j++)
 512                 tmp_coi[j] = *((__global short*)src_data+(y+lidy)*src_step/2+(x+i+j)*kcn+coi-1);
 513             tmp[i/VLEN_S] = (short8)(tmp_coi[0],tmp_coi[1],tmp_coi[2],tmp_coi[3],tmp_coi[4],tmp_coi[5],tmp_coi[6],tmp_coi[7]);
 514         }
 515     else
 516         for(int i=0; i < tileSize_width; i+=VLEN_S)
 517             tmp[i/VLEN_S] = *(src_data+(y+lidy)*src_step/(2*VLEN_S)+(x+i)/VLEN_S);
 518     short8 zero = (short8)(0);
 519     short8 full = (short8)(255);
 520     if( binary )
 521         for(int i=0; i < tileSize_width; i+=(VLEN_S))
 522             tmp[i/VLEN_S] = (tmp[i/VLEN_S]!=zero)?full:zero;
 523
 524     F mom[10];
 525     __local long m[10][128];
 526     if(lidy == 0)
 527         for(int i=0; i<10; i++)
 528             for(int j=0; j<128; j++)
 529                 m[i][j]=0;
 530     barrier(CLK_LOCAL_MEM_FENCE);
 531     long lm[10] = {0};
 532     int8 x0 = (int8)(0);
 533     int8 x1 = (int8)(0);
 534     int8 x2 = (int8)(0);
 535     long8 x3 = (long8)(0);
 536     for( int xt = 0 ; xt < tileSize_width; xt+= (VLEN_S))
 537     {
 538         int8 v_xt = (int8)(xt, xt+1, xt+2, xt+3, xt+4, xt+5, xt+6, xt+7);
 539         int8 p = convert_int8(tmp[xt/VLEN_S]);
 540         int8 xp = v_xt * p, xxp = xp * v_xt;
 541         x0 += p;
 542         x1 += xp;
 543         x2 += xxp;
 544         x3 += convert_long8(xxp) * convert_long8(v_xt);
 545     }
 546     x0.s0 += x0.s1 + x0.s2 + x0.s3 + x0.s4 + x0.s5 + x0.s6 + x0.s7;
 547     x1.s0 += x1.s1 + x1.s2 + x1.s3 + x1.s4 + x1.s5 + x1.s6 + x1.s7;
 548     x2.s0 += x2.s1 + x2.s2 + x2.s3 + x2.s4 + x2.s5 + x2.s6 + x2.s7;
 549     x3.s0 += x3.s1 + x3.s2 + x3.s3 + x3.s4 + x3.s5 + x3.s6 + x3.s7;
 550
 551     int py = lidy * x0.s0, sy = lidy*lidy;
 552     int bheight = min(tileSize_height, TILE_SIZE/2);
 553     if(bheight >= TILE_SIZE/2&&lidy > bheight-1&&lidy < tileSize_height)
 554     {
 555         m[9][lidy-bheight] = ((long)py) * sy;  // m03
 556         m[8][lidy-bheight] = ((long)x1.s0) * sy;  // m12
 557         m[7][lidy-bheight] = ((long)x2.s0) * lidy;  // m21
 558         m[6][lidy-bheight] = x3.s0;             // m30
 559         m[5][lidy-bheight] = x0.s0 * sy;        // m02
 560         m[4][lidy-bheight] = x1.s0 * lidy;         // m11
 561         m[3][lidy-bheight] = x2.s0;             // m20
 562         m[2][lidy-bheight] = py;             // m01
 563         m[1][lidy-bheight] = x1.s0;             // m10
 564         m[0][lidy-bheight] = x0.s0;             // m00
 565     }
 566     else if(lidy < bheight)
 567     {
 568         lm[9] = ((long)py) * sy;  // m03
 569         lm[8] = ((long)(x1.s0)) * sy;  // m12
 570         lm[7] = ((long)(x2.s0)) * lidy;  // m21
 571         lm[6] = x3.s0;             // m30
 572         lm[5] = x0.s0 * sy;        // m02
 573         lm[4] = x1.s0 * lidy;         // m11
 574         lm[3] = x2.s0;             // m20
 575         lm[2] = py;             // m01
 576         lm[1] = x1.s0;             // m10
 577         lm[0] = x0.s0;             // m00
 578     }
 579     barrier(CLK_LOCAL_MEM_FENCE);
 580     for( int j = TILE_SIZE/2; j >=1; j = j/2 )
 581     {
 582         if(lidy < j)
 583             for( int i = 0; i < 10; i++ )
 584                 lm[i] = lm[i] + m[i][lidy];
 585         barrier(CLK_LOCAL_MEM_FENCE);
 586         if(lidy >= j/2&&lidy < j)
 587             for( int i = 0; i < 10; i++ )
 588                 m[i][lidy-j/2] = lm[i];
 589         barrier(CLK_LOCAL_MEM_FENCE);
 590     }
 591     if(lidy ==0 &&lidx ==0)
 592     {
 593         for(int mt = 0; mt < 10; mt++ )
 594             mom[mt] = (F)lm[mt];
 595
 596         if(binary)
 597         {
 598             F s = 1./255;
 599             for( int mt = 0; mt < 10; mt++ )
 600                 mom[mt] *= s;
 601         }
 602
 603         F xm = x * mom[0], ym = y*mom[0];
 604
 605         // accumulate moments computed in each tile
 606         dst_step /= sizeof(F);
 607
 608         // + m00 ( = m00' )
 609         *(dst_m + mad24(DST_ROW_00 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[0];
 610
 611         // + m10 ( = m10' + x*m00' )
 612         *(dst_m + mad24(DST_ROW_10 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[1] + xm;
 613
 614         // + m01 ( = m01' + y*m00' )
 615         *(dst_m + mad24(DST_ROW_01 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[2] + ym;
 616
 617         // + m20 ( = m20' + 2*x*m10' + x*x*m00' )
 618         *(dst_m + mad24(DST_ROW_20 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[3] + x * (mom[1] * 2 + xm);
 619
 620         // + m11 ( = m11' + x*m01' + y*m10' + x*y*m00' )
 621         *(dst_m + mad24(DST_ROW_11 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[4] + x * (mom[2] + ym) + y * mom[1];
 622
 623         // + m02 ( = m02' + 2*y*m01' + y*y*m00' )
 624         *(dst_m + mad24(DST_ROW_02 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[5] + y * (mom[2] * 2 + ym);
 625
 626         // + m30 ( = m30' + 3*x*m20' + 3*x*x*m10' + x*x*x*m00' )
 627         *(dst_m + mad24(DST_ROW_30 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[6] + x * (3. * mom[3] + x * (3. * mom[1] + xm));
 628
 629         // + m21 ( = m21' + x*(2*m11' + 2*y*m10' + x*m01' + x*y*m00') + y*m20')
 630         *(dst_m + mad24(DST_ROW_21 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[7] + x * (2 * (mom[4] + y * mom[1]) + x * (mom[2] + ym)) + y * mom[3];
 631
 632         // + m12 ( = m12' + y*(2*m11' + 2*x*m01' + y*m10' + x*y*m00') + x*m02')
 633         *(dst_m + mad24(DST_ROW_12 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[8] + y * (2 * (mom[4] + x * mom[2]) + y * (mom[1] + xm)) + x * mom[5];
 634
 635         // + m03 ( = m03' + 3*y*m02' + 3*y*y*m01' + y*y*y*m00' )
 636         *(dst_m + mad24(DST_ROW_03 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[9] + y * (3. * mom[5] + y * (3. * mom[2] + ym));
 637     }
 638 }
 639
 640 __kernel void CvMoments_D5( __global float* src_data, int src_rows, int src_cols, int src_step, int tileSize_width, int tileSize_height,
 641                             __global F* dst_m,
 642                             int dst_cols, int dst_step, int blocky,
 643                             int type, int depth, int cn, int coi, int binary, const int TILE_SIZE)
 644 {
 645     float tmp_coi[4]; // get the coi data
 646     float4 tmp[64] ;
 647     int VLEN_F = 4; // vector length of float
 648     int gidy = get_global_id(0);
 649     int gidx = get_global_id(1);
 650     int wgidy = get_group_id(0);
 651     int wgidx = get_group_id(1);
 652     int lidy = get_local_id(0);
 653     int lidx = get_local_id(1);
 654     int y = wgidy*TILE_SIZE;  // real Y index of pixel
 655     int x = wgidx*TILE_SIZE;  // real X index of pixel
 656     int kcn = (cn==2)?2:4;
 657     src_step /= sizeof(*src_data);
 658     int rstep = min(src_step, TILE_SIZE);
 659     tileSize_height = min(TILE_SIZE, src_rows - y);
 660     tileSize_width = min(TILE_SIZE, src_cols -x);
 661     int maxIdx = mul24(src_rows, src_cols);
 662     int yOff = (y+lidy)*src_step;
 663     int index;
 664     if(tileSize_width < TILE_SIZE && yOff < src_rows)
 665         for(int i = tileSize_width; i < rstep && (yOff+x+i) < maxIdx; i++ )
 666             *(src_data+yOff+x+i) = 0;
 667     if( coi > 0 )
 668         for(int i=0; i < tileSize_width; i+=VLEN_F)
 669         {
 670 #pragma unroll
 671             for(int j=0; j<4; j++)
 672             {
 673                 index = yOff+(x+i+j)*kcn+coi-1;
 674                 if (index < maxIdx)
 675                     tmp_coi[j] = *(src_data+index);
 676                 else
 677                     tmp_coi[j] = 0;
 678             }
 679             tmp[i/VLEN_F] = (float4)(tmp_coi[0],tmp_coi[1],tmp_coi[2],tmp_coi[3]);
 680         }
 681     else
 682         for(int i=0; i < tileSize_width && (yOff+x+i) < maxIdx; i+=VLEN_F)
 683             tmp[i/VLEN_F] = (*(__global float4 *)(src_data+yOff+x+i));
 684     float4 zero = (float4)(0);
 685     float4 full = (float4)(255);
 686     if( binary )
 687         for(int i=0; i < tileSize_width; i+=4)
 688             tmp[i/VLEN_F] = (tmp[i/VLEN_F]!=zero)?full:zero;
 689     F mom[10];
 690     __local F m[10][128];
 691     if(lidy == 0)
 692         for(int i = 0; i < 10; i ++)
 693             for(int j = 0; j < 128; j ++)
 694                 m[i][j] = 0;
 695     barrier(CLK_LOCAL_MEM_FENCE);
 696     F lm[10] = {0};
 697     F4 x0 = (F4)(0);
 698     F4 x1 = (F4)(0);
 699     F4 x2 = (F4)(0);
 700     F4 x3 = (F4)(0);
 701     for( int xt = 0 ; xt < tileSize_width; xt+=VLEN_F )
 702     {
 703         F4 v_xt = (F4)(xt, xt+1, xt+2, xt+3);
 704         F4 p = convert_F4(tmp[xt/VLEN_F]);
 705         F4 xp = v_xt * p, xxp = xp * v_xt;
 706         x0 += p;
 707         x1 += xp;
 708         x2 += xxp;
 709         x3 += xxp * v_xt;
 710     }
 711     x0.s0 += x0.s1 + x0.s2 + x0.s3;
 712     x1.s0 += x1.s1 + x1.s2 + x1.s3;
 713     x2.s0 += x2.s1 + x2.s2 + x2.s3;
 714     x3.s0 += x3.s1 + x3.s2 + x3.s3;
 715
 716     F py = lidy * x0.s0, sy = lidy*lidy;
 717     int bheight = min(tileSize_height, TILE_SIZE/2);
 718     if(bheight >= TILE_SIZE/2&&lidy > bheight-1&&lidy < tileSize_height)
 719     {
 720         m[9][lidy-bheight] = ((F)py) * sy;  // m03
 721         m[8][lidy-bheight] = ((F)x1.s0) * sy;  // m12
 722         m[7][lidy-bheight] = ((F)x2.s0) * lidy;  // m21
 723         m[6][lidy-bheight] = x3.s0;             // m30
 724         m[5][lidy-bheight] = x0.s0 * sy;        // m02
 725         m[4][lidy-bheight] = x1.s0 * lidy;         // m11
 726         m[3][lidy-bheight] = x2.s0;             // m20
 727         m[2][lidy-bheight] = py;             // m01
 728         m[1][lidy-bheight] = x1.s0;             // m10
 729         m[0][lidy-bheight] = x0.s0;             // m00
 730     }
 731
 732     else if(lidy < bheight)
 733     {
 734         lm[9] = ((F)py) * sy;  // m03
 735         lm[8] = ((F)x1.s0) * sy;  // m12
 736         lm[7] = ((F)x2.s0) * lidy;  // m21
 737         lm[6] = x3.s0;             // m30
 738         lm[5] = x0.s0 * sy;        // m02
 739         lm[4] = x1.s0 * lidy;         // m11
 740         lm[3] = x2.s0;             // m20
 741         lm[2] = py;             // m01
 742         lm[1] = x1.s0;             // m10
 743         lm[0] = x0.s0;             // m00
 744     }
 745     barrier(CLK_LOCAL_MEM_FENCE);
 746     for( int j = TILE_SIZE/2; j >= 1; j = j/2 )
 747     {
 748         if(lidy < j)
 749             for( int i = 0; i < 10; i++ )
 750                 lm[i] = lm[i] + m[i][lidy];
 751         barrier(CLK_LOCAL_MEM_FENCE);
 752         if(lidy >= j/2&&lidy < j)
 753             for( int i = 0; i < 10; i++ )
 754                 m[i][lidy-j/2] = lm[i];
 755         barrier(CLK_LOCAL_MEM_FENCE);
 756     }
 757     if(lidy == 0&&lidx == 0)
 758     {
 759         for( int mt = 0; mt < 10; mt++ )
 760             mom[mt] = (F)lm[mt];
 761         if(binary)
 762         {
 763             F s = 1./255;
 764             for( int mt = 0; mt < 10; mt++ )
 765                 mom[mt] *= s;
 766         }
 767
 768         F xm = x * mom[0], ym = y * mom[0];
 769
 770         // accumulate moments computed in each tile
 771         dst_step /= sizeof(F);
 772
 773         int dst_x_off = mad24(wgidy, dst_cols, wgidx);
 774         int dst_off = 0;
 775         int max_dst_index = 10 * blocky * get_global_size(1);
 776
 777         // + m00 ( = m00' )
 778         dst_off = mad24(DST_ROW_00 * blocky, dst_step, dst_x_off);
 779         if (dst_off < max_dst_index)
 780             *(dst_m + dst_off) = mom[0];
 781
 782         // + m10 ( = m10' + x*m00' )
 783         dst_off = mad24(DST_ROW_10 * blocky, dst_step, dst_x_off);
 784         if (dst_off < max_dst_index)
 785             *(dst_m + dst_off) = mom[1] + xm;
 786
 787         // + m01 ( = m01' + y*m00' )
 788         dst_off = mad24(DST_ROW_01 * blocky, dst_step, dst_x_off);
 789         if (dst_off < max_dst_index)
 790             *(dst_m + dst_off) = mom[2] + ym;
 791
 792         // + m20 ( = m20' + 2*x*m10' + x*x*m00' )
 793         dst_off = mad24(DST_ROW_20 * blocky, dst_step, dst_x_off);
 794         if (dst_off < max_dst_index)
 795             *(dst_m + dst_off) = mom[3] + x * (mom[1] * 2 + xm);
 796
 797         // + m11 ( = m11' + x*m01' + y*m10' + x*y*m00' )
 798         dst_off = mad24(DST_ROW_11 * blocky, dst_step, dst_x_off);
 799         if (dst_off < max_dst_index)
 800             *(dst_m + dst_off) = mom[4] + x * (mom[2] + ym) + y * mom[1];
 801
 802         // + m02 ( = m02' + 2*y*m01' + y*y*m00' )
 803         dst_off = mad24(DST_ROW_02 * blocky, dst_step, dst_x_off);
 804         if (dst_off < max_dst_index)
 805             *(dst_m + dst_off) = mom[5] + y * (mom[2] * 2 + ym);
 806
 807         // + m30 ( = m30' + 3*x*m20' + 3*x*x*m10' + x*x*x*m00' )
 808         dst_off = mad24(DST_ROW_30 * blocky, dst_step, dst_x_off);
 809         if (dst_off < max_dst_index)
 810             *(dst_m + dst_off) = mom[6] + x * (3. * mom[3] + x * (3. * mom[1] + xm));
 811
 812         // + m21 ( = m21' + x*(2*m11' + 2*y*m10' + x*m01' + x*y*m00') + y*m20')
 813         dst_off = mad24(DST_ROW_21 * blocky, dst_step, dst_x_off);
 814         if (dst_off < max_dst_index)
 815             *(dst_m + dst_off) = mom[7] + x * (2 * (mom[4] + y * mom[1]) + x * (mom[2] + ym)) + y * mom[3];
 816
 817         // + m12 ( = m12' + y*(2*m11' + 2*x*m01' + y*m10' + x*y*m00') + x*m02')
 818         dst_off = mad24(DST_ROW_12 * blocky, dst_step, dst_x_off);
 819         if (dst_off < max_dst_index)
 820             *(dst_m + dst_off) = mom[8] + y * (2 * (mom[4] + x * mom[2]) + y * (mom[1] + xm)) + x * mom[5];
 821
 822         // + m03 ( = m03' + 3*y*m02' + 3*y*y*m01' + y*y*y*m00' )
 823         dst_off = mad24(DST_ROW_03 * blocky, dst_step, dst_x_off);
 824         if (dst_off < max_dst_index)
 825             *(dst_m + dst_off) = mom[9] + y * (3. * mom[5] + y * (3. * mom[2] + ym));
 826     }
 827 }
 828
 829 __kernel void CvMoments_D6(__global F* src_data,  int src_rows, int src_cols, int src_step, int tileSize_width, int tileSize_height,
 830                            __global F* dst_m,
 831                            int dst_cols, int dst_step, int blocky,
 832                            int type, int depth, int cn, int coi, int binary, const int TILE_SIZE)
 833 {
 834     F tmp_coi[4]; // get the coi data
 835     F4 tmp[64];
 836     int VLEN_D = 4; // length of vetor
 837     int gidy = get_global_id(0);
 838     int gidx = get_global_id(1);
 839     int wgidy = get_group_id(0);
 840     int wgidx = get_group_id(1);
 841     int lidy = get_local_id(0);
 842     int lidx = get_local_id(1);
 843     int y = wgidy*TILE_SIZE;  // real Y index of pixel
 844     int x = wgidx*TILE_SIZE;  // real X index of pixel
 845     int kcn = (cn==2)?2:4;
 846     int rstep = min(src_step/8, TILE_SIZE);
 847     tileSize_height = min(TILE_SIZE,  src_rows - y);
 848     tileSize_width = min(TILE_SIZE, src_cols - x);
 849
 850     if(tileSize_width < TILE_SIZE)
 851         for(int i = tileSize_width; i < rstep; i++ )
 852             *((__global F*)src_data+(y+lidy)*src_step/8+x+i) = 0;
 853     if( coi > 0 )
 854         for(int i=0; i < tileSize_width; i+=VLEN_D)
 855         {
 856             for(int j=0; j<4; j++)
 857                 tmp_coi[j] = *(src_data+(y+lidy)*src_step/8+(x+i+j)*kcn+coi-1);
 858             tmp[i/VLEN_D] = (F4)(tmp_coi[0],tmp_coi[1],tmp_coi[2],tmp_coi[3]);
 859         }
 860     else
 861         for(int i=0; i < tileSize_width; i+=VLEN_D)
 862             tmp[i/VLEN_D] = (F4)(*(src_data+(y+lidy)*src_step/8+x+i),*(src_data+(y+lidy)*src_step/8+x+i+1),*(src_data+(y+lidy)*src_step/8+x+i+2),*(src_data+(y+lidy)*src_step/8+x+i+3));
 863     F4 zero = (F4)(0);
 864     F4 full = (F4)(255);
 865     if( binary )
 866         for(int i=0; i < tileSize_width; i+=VLEN_D)
 867             tmp[i/VLEN_D] = (tmp[i/VLEN_D]!=zero)?full:zero;
 868     F mom[10];
 869     __local F m[10][128];
 870     if(lidy == 0)
 871         for(int i=0; i<10; i++)
 872             for(int j=0; j<128; j++)
 873                 m[i][j]=0;
 874     barrier(CLK_LOCAL_MEM_FENCE);
 875     F lm[10] = {0};
 876     F4 x0 = (F4)(0);
 877     F4 x1 = (F4)(0);
 878     F4 x2 = (F4)(0);
 879     F4 x3 = (F4)(0);
 880     for( int xt = 0 ; xt < tileSize_width; xt+=VLEN_D )
 881     {
 882         F4 v_xt = (F4)(xt, xt+1, xt+2, xt+3);
 883         F4 p = tmp[xt/VLEN_D];
 884         F4 xp = v_xt * p, xxp = xp * v_xt;
 885         x0 += p;
 886         x1 += xp;
 887         x2 += xxp;
 888         x3 += xxp *v_xt;
 889     }
 890     x0.s0 += x0.s1 + x0.s2 + x0.s3;
 891     x1.s0 += x1.s1 + x1.s2 + x1.s3;
 892     x2.s0 += x2.s1 + x2.s2 + x2.s3;
 893     x3.s0 += x3.s1 + x3.s2 + x3.s3;
 894
 895     F py = lidy * x0.s0, sy = lidy*lidy;
 896     int bheight = min(tileSize_height, TILE_SIZE/2);
 897     if(bheight >= TILE_SIZE/2&&lidy > bheight-1&&lidy < tileSize_height)
 898     {
 899         m[9][lidy-bheight] = ((F)py) * sy;  // m03
 900         m[8][lidy-bheight] = ((F)x1.s0) * sy;  // m12
 901         m[7][lidy-bheight] = ((F)x2.s0) * lidy;  // m21
 902         m[6][lidy-bheight] = x3.s0;             // m30
 903         m[5][lidy-bheight] = x0.s0 * sy;        // m02
 904         m[4][lidy-bheight] = x1.s0 * lidy;         // m11
 905         m[3][lidy-bheight] = x2.s0;             // m20
 906         m[2][lidy-bheight] = py;             // m01
 907         m[1][lidy-bheight] = x1.s0;             // m10
 908         m[0][lidy-bheight] = x0.s0;             // m00
 909     }
 910
 911     else if(lidy < bheight)
 912     {
 913         lm[9] = ((F)py) * sy;  // m03
 914         lm[8] = ((F)x1.s0) * sy;  // m12
 915         lm[7] = ((F)x2.s0) * lidy;  // m21
 916         lm[6] = x3.s0;             // m30
 917         lm[5] = x0.s0 * sy;        // m02
 918         lm[4] = x1.s0 * lidy;         // m11
 919         lm[3] = x2.s0;             // m20
 920         lm[2] = py;             // m01
 921         lm[1] = x1.s0;             // m10
 922         lm[0] = x0.s0;             // m00
 923     }
 924     barrier(CLK_LOCAL_MEM_FENCE);
 925     for( int j = TILE_SIZE/2; j >= 1; j = j/2 )
 926     {
 927         if(lidy < j)
 928             for( int i = 0; i < 10; i++ )
 929                 lm[i] = lm[i] + m[i][lidy];
 930         barrier(CLK_LOCAL_MEM_FENCE);
 931         if(lidy >= j/2&&lidy < j)
 932             for( int i = 0; i < 10; i++ )
 933                 m[i][lidy-j/2] = lm[i];
 934         barrier(CLK_LOCAL_MEM_FENCE);
 935     }
 936     if(lidy == 0&&lidx == 0)
 937     {
 938         for( int mt = 0; mt < 10; mt++ )
 939             mom[mt] = (F)lm[mt];
 940         if(binary)
 941         {
 942             F s = 1./255;
 943             for( int mt = 0; mt < 10; mt++ )
 944                 mom[mt] *= s;
 945         }
 946
 947         F xm = x * mom[0], ym = y * mom[0];
 948
 949         // accumulate moments computed in each tile
 950         dst_step /= sizeof(F);
 951
 952         // + m00 ( = m00' )
 953         *(dst_m + mad24(DST_ROW_00 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[0];
 954
 955         // + m10 ( = m10' + x*m00' )
 956         *(dst_m + mad24(DST_ROW_10 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[1] + xm;
 957
 958         // + m01 ( = m01' + y*m00' )
 959         *(dst_m + mad24(DST_ROW_01 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[2] + ym;
 960
 961         // + m20 ( = m20' + 2*x*m10' + x*x*m00' )
 962         *(dst_m + mad24(DST_ROW_20 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[3] + x * (mom[1] * 2 + xm);
 963
 964         // + m11 ( = m11' + x*m01' + y*m10' + x*y*m00' )
 965         *(dst_m + mad24(DST_ROW_11 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[4] + x * (mom[2] + ym) + y * mom[1];
 966
 967         // + m02 ( = m02' + 2*y*m01' + y*y*m00' )
 968         *(dst_m + mad24(DST_ROW_02 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[5] + y * (mom[2] * 2 + ym);
 969
 970         // + m30 ( = m30' + 3*x*m20' + 3*x*x*m10' + x*x*x*m00' )
 971         *(dst_m + mad24(DST_ROW_30 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[6] + x * (3. * mom[3] + x * (3. * mom[1] + xm));
 972
 973         // + m21 ( = m21' + x*(2*m11' + 2*y*m10' + x*m01' + x*y*m00') + y*m20')
 974         *(dst_m + mad24(DST_ROW_21 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[7] + x * (2 * (mom[4] + y * mom[1]) + x * (mom[2] + ym)) + y * mom[3];
 975
 976         // + m12 ( = m12' + y*(2*m11' + 2*x*m01' + y*m10' + x*y*m00') + x*m02')
 977         *(dst_m + mad24(DST_ROW_12 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[8] + y * (2 * (mom[4] + x * mom[2]) + y * (mom[1] + xm)) + x * mom[5];
 978
 979         // + m03 ( = m03' + 3*y*m02' + 3*y*y*m01' + y*y*y*m00' )
 980         *(dst_m + mad24(DST_ROW_03 * blocky, dst_step, mad24(wgidy, dst_cols, wgidx))) = mom[9] + y * (3. * mom[5] + y * (3. * mom[2] + ym));
 981     }
 982 }