A few optimizations to ocl::pyrLK::sparse, make it running on more OCL platforms
authoryao <bitwangyaoyao@gmail.com>
Wed, 16 Jan 2013 08:50:59 +0000 (16:50 +0800)
committerAndrey Kamaev <andrey.kamaev@itseez.com>
Wed, 23 Jan 2013 10:48:04 +0000 (14:48 +0400)
modules/ocl/src/kernels/arithm_mul.cl
modules/ocl/src/kernels/pyrlk.cl
modules/ocl/src/kernels/pyrlk_no_image.cl [new file with mode: 0644]
modules/ocl/src/pyrlk.cpp
modules/ocl/test/test_pyrlk.cpp

index e0cfbd8..f9f3936 100644 (file)
@@ -16,6 +16,7 @@
 //
 // @Authors
 //    Jia Haipeng, jiahaipeng95@gmail.com
+//    Dachuan Zhao, dachuan@multicorewareinc.com
 //
 // Redistribution and use in source and binary forms, with or without modification,
 // are permitted provided that the following conditions are met:
@@ -260,3 +261,22 @@ __kernel void arithm_mul_D6 (__global double *src1, int src1_step, int src1_offs
     }
 }
 #endif
+
+__kernel void arithm_muls_D5 (__global float *src1, int src1_step, int src1_offset,
+                             __global float *dst,  int dst_step,  int dst_offset,
+                             int rows, int cols, int dst_step1, float scalar)
+{
+    int x = get_global_id(0);
+    int y = get_global_id(1);
+
+    if (x < cols && y < rows)
+    {
+        int src1_index = mad24(y, src1_step, (x << 2) + src1_offset);
+        int dst_index  = mad24(y, dst_step,  (x << 2) + dst_offset);
+
+        float data1 = *((__global float *)((__global char *)src1 + src1_index));
+        float tmp = data1 * scalar;
+
+        *((__global float *)((__global char *)dst + dst_index)) = tmp;
+    }
+}
\ No newline at end of file
index ecdacf3..c772be7 100644 (file)
@@ -16,6 +16,7 @@
 //
 // @Authors
 //    Dachuan Zhao, dachuan@multicorewareinc.com
+//    Yao Wang, bitwangyaoyao@gmail.com
 //
 // Redistribution and use in source and binary forms, with or without modification,
 // are permitted provided that the following conditions are met:
 
 //#pragma OPENCL EXTENSION cl_amd_printf : enable
 
-__kernel void arithm_muls_D5 (__global float *src1, int src1_step, int src1_offset,
-                             __global float *dst,  int dst_step,  int dst_offset,
-                             int rows, int cols, int dst_step1, float scalar)
-{
-    int x = get_global_id(0);
-    int y = get_global_id(1);
-
-    if (x < cols && y < rows)
-    {
-        int src1_index = mad24(y, src1_step, (x << 2) + src1_offset);
-        int dst_index  = mad24(y, dst_step,  (x << 2) + dst_offset);
-
-        float data1 = *((__global float *)((__global char *)src1 + src1_index));
-        float tmp = data1 * scalar;
-
-        *((__global float *)((__global char *)dst + dst_index)) = tmp;
-    }
-}
-
-
 __kernel void calcSharrDeriv_vertical_C1_D0(__global const uchar* src, int srcStep, int rows, int cols, int cn, __global short* dx_buf, int dx_bufStep, __global short* dy_buf, int dy_bufStep)
 {
     const int x = get_global_id(0);
@@ -202,6 +183,7 @@ float linearFilter_float(__global const float* src, int srcStep, int cn, float2
     return src_row[(int)x] * iw00 + src_row[(int)x + cn] * iw01 + src_row1[(int)x] * iw10 + src_row1[(int)x + cn] * iw11, W_BITS1 - 5;
 }
 
+#define        BUFFER  64
 void reduce3(float val1, float val2, float val3, __local float* smem1, __local float* smem2, __local float* smem3, int tid)
 {
     smem1[tid] = val1;
@@ -209,6 +191,7 @@ void reduce3(float val1, float val2, float val3, __local float* smem1, __local f
     smem3[tid] = val3;
     barrier(CLK_LOCAL_MEM_FENCE);
 
+#if    BUFFER > 128
     if (tid < 128)
     {
         smem1[tid] = val1 += smem1[tid + 128];
@@ -216,7 +199,9 @@ void reduce3(float val1, float val2, float val3, __local float* smem1, __local f
         smem3[tid] = val3 += smem3[tid + 128];
     }
     barrier(CLK_LOCAL_MEM_FENCE);
+#endif
 
+#if    BUFFER > 64
     if (tid < 64)
     {
         smem1[tid] = val1 += smem1[tid + 64];
@@ -224,6 +209,7 @@ void reduce3(float val1, float val2, float val3, __local float* smem1, __local f
         smem3[tid] = val3 += smem3[tid + 64];
     }
     barrier(CLK_LOCAL_MEM_FENCE);
+#endif
 
     if (tid < 32)
     {
@@ -263,19 +249,23 @@ void reduce2(float val1, float val2, __local float* smem1, __local float* smem2,
     smem2[tid] = val2;
     barrier(CLK_LOCAL_MEM_FENCE);
 
+#if    BUFFER > 128
     if (tid < 128)
     {
         smem1[tid] = val1 += smem1[tid + 128];
         smem2[tid] = val2 += smem2[tid + 128];
     }
     barrier(CLK_LOCAL_MEM_FENCE);
+#endif
 
+#if    BUFFER > 64
     if (tid < 64)
     {
         smem1[tid] = val1 += smem1[tid + 64];
         smem2[tid] = val2 += smem2[tid + 64];
     }
     barrier(CLK_LOCAL_MEM_FENCE);
+#endif
 
     if (tid < 32)
     {
@@ -307,17 +297,21 @@ void reduce1(float val1, __local float* smem1, int tid)
     smem1[tid] = val1;
     barrier(CLK_LOCAL_MEM_FENCE);
 
+#if    BUFFER > 128
     if (tid < 128)
     {
         smem1[tid] = val1 += smem1[tid + 128];
     }
     barrier(CLK_LOCAL_MEM_FENCE);
+#endif
 
+#if    BUFFER > 64
     if (tid < 64)
     {
         smem1[tid] = val1 += smem1[tid + 64];
     }
     barrier(CLK_LOCAL_MEM_FENCE);
+#endif
 
     if (tid < 32)
     {
@@ -333,41 +327,122 @@ void reduce1(float val1, __local float* smem1, int tid)
 }
 
 #define SCALE (1.0f / (1 << 20))
+#define        THRESHOLD       0.01f
+#define        DIMENSION       21
 
 // Image read mode
 __constant sampler_t sampler    = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_LINEAR;
 
+void SetPatch(image2d_t I, float x, float y,
+                                float* Pch, float* Dx, float* Dy,
+                                float* A11, float* A12, float* A22)
+{
+            *Pch = read_imagef(I, sampler, (float2)(x, y)).x;
+
+            float dIdx = 3.0f * read_imagef(I, sampler, (float2)(x + 1, y - 1)).x + 10.0f * read_imagef(I, sampler, (float2)(x + 1, y)).x + 3.0f * read_imagef(I, sampler, (float2)(x + 1, y + 1)).x -
+                             (3.0f * read_imagef(I, sampler, (float2)(x - 1, y - 1)).x + 10.0f * read_imagef(I, sampler, (float2)(x - 1, y)).x + 3.0f * read_imagef(I, sampler, (float2)(x - 1, y + 1)).x);
+
+            float dIdy = 3.0f * read_imagef(I, sampler, (float2)(x - 1, y + 1)).x + 10.0f * read_imagef(I, sampler, (float2)(x, y + 1)).x + 3.0f * read_imagef(I, sampler, (float2)(x + 1, y + 1)).x -
+                            (3.0f * read_imagef(I, sampler, (float2)(x - 1, y - 1)).x + 10.0f * read_imagef(I, sampler, (float2)(x, y - 1)).x + 3.0f * read_imagef(I, sampler, (float2)(x + 1, y - 1)).x);
+
+
+            *Dx = dIdx;
+            *Dy = dIdy;
+
+            *A11 += dIdx * dIdx;
+            *A12 += dIdx * dIdy;
+            *A22 += dIdy * dIdy;
+}
+
+void GetPatch(image2d_t J, float x, float y,
+                                float* Pch, float* Dx, float* Dy,
+                                float* b1, float* b2)
+{
+                float J_val = read_imagef(J, sampler, (float2)(x, y)).x;
+                float diff = (J_val - *Pch) * 32.0f;
+                *b1 += diff**Dx;
+                *b2 += diff**Dy;
+}
+
+void GetError(image2d_t J, const float x, const float y, const float* Pch, float* errval)
+{
+        float diff = read_imagef(J, sampler, (float2)(x,y)).x-*Pch;
+        *errval += fabs(diff);
+}
+
+void SetPatch4(image2d_t I, const float x, const float y,
+                                float4* Pch, float4* Dx, float4* Dy,
+                                float* A11, float* A12, float* A22)
+{
+            *Pch = read_imagef(I, sampler, (float2)(x, y));
+
+            float4 dIdx = 3.0f * read_imagef(I, sampler, (float2)(x + 1, y - 1)) + 10.0f * read_imagef(I, sampler, (float2)(x + 1, y)) + 3.0f * read_imagef(I, sampler, (float2)(x + 1, y + 1)) -
+                             (3.0f * read_imagef(I, sampler, (float2)(x - 1, y - 1)) + 10.0f * read_imagef(I, sampler, (float2)(x - 1, y)) + 3.0f * read_imagef(I, sampler, (float2)(x - 1, y + 1)));
+
+            float4 dIdy = 3.0f * read_imagef(I, sampler, (float2)(x - 1, y + 1)) + 10.0f * read_imagef(I, sampler, (float2)(x, y + 1)) + 3.0f * read_imagef(I, sampler, (float2)(x + 1, y + 1)) -
+                            (3.0f * read_imagef(I, sampler, (float2)(x - 1, y - 1)) + 10.0f * read_imagef(I, sampler, (float2)(x, y - 1)) + 3.0f * read_imagef(I, sampler, (float2)(x + 1, y - 1)));
+
+
+            *Dx = dIdx;
+            *Dy = dIdy;
+                        float4 sqIdx = dIdx * dIdx;
+                        *A11 += sqIdx.x + sqIdx.y + sqIdx.z;
+                        sqIdx = dIdx * dIdy;
+                        *A12 += sqIdx.x + sqIdx.y + sqIdx.z;
+                        sqIdx = dIdy * dIdy;
+                        *A22 += sqIdx.x + sqIdx.y + sqIdx.z;
+}
+
+void GetPatch4(image2d_t J, const float x, const float y,
+                                const float4* Pch, const float4* Dx, const float4* Dy,
+                                float* b1, float* b2)
+{
+                float4 J_val = read_imagef(J, sampler, (float2)(x, y));
+                float4 diff = (J_val - *Pch) * 32.0f;
+                                float4 xdiff = diff* *Dx;
+                                *b1 += xdiff.x + xdiff.y + xdiff.z;
+                                xdiff = diff* *Dy;
+                                *b2 += xdiff.x + xdiff.y + xdiff.z;
+}
+
+void GetError4(image2d_t J, const float x, const float y, const float4* Pch, float* errval)
+{
+        float4 diff = read_imagef(J, sampler, (float2)(x,y))-*Pch;
+        *errval += fabs(diff.x) + fabs(diff.y) + fabs(diff.z);
+}
+
+
 __kernel void lkSparse_C1_D5(image2d_t I, image2d_t J,
-    __global const float2* prevPts, int prevPtsStep, __global float2* nextPts, int nextPtsStep, __global uchar* status/*, __global float* err*/, const int level, const int rows, const int cols, int PATCH_X, int PATCH_Y, int cn, int c_winSize_x, int c_winSize_y, int c_iters, char calcErr, char GET_MIN_EIGENVALS)
+    __global const float2* prevPts, int prevPtsStep, __global float2* nextPts, int nextPtsStep, __global uchar* status, __global float* err,
+        const int level, const int rows, const int cols, int PATCH_X, int PATCH_Y, int cn, int c_winSize_x, int c_winSize_y, int c_iters, char calcErr)
 {
-    __local float smem1[256];
-    __local float smem2[256];
-    __local float smem3[256];
+    __local float smem1[BUFFER];
+    __local float smem2[BUFFER];
+    __local float smem3[BUFFER];
+
+        unsigned int xid=get_local_id(0);
+        unsigned int yid=get_local_id(1);
+        unsigned int gid=get_group_id(0);
+        unsigned int xsize=get_local_size(0);
+        unsigned int ysize=get_local_size(1);
+        int xBase, yBase, i, j, k;
 
-    int c_halfWin_x = (c_winSize_x - 1) / 2;
-    int c_halfWin_y = (c_winSize_y - 1) / 2;
+        float2 c_halfWin = (float2)((c_winSize_x - 1)>>1, (c_winSize_y - 1)>>1);
 
-    const int tid = get_local_id(1) * get_local_size(0) + get_local_id(0);
+    const int tid = mad24(yid, xsize, xid);
 
-    float2 prevPt = prevPts[get_group_id(0)];
-    prevPt.x *= (1.0f / (1 << level));
-    prevPt.y *= (1.0f / (1 << level));
+    float2 prevPt = prevPts[gid] / (1 << level);
 
     if (prevPt.x < 0 || prevPt.x >= cols || prevPt.y < 0 || prevPt.y >= rows)
     {
-        if (level == 0 && tid == 0)
+        if (tid == 0 && level == 0)
         {
-            status[get_group_id(0)] = 0;
-
-            //if (calcErr)
-            //    err[get_group_id(0)] = 0;
+            status[gid] = 0;
         }
 
         return;
     }
-
-    prevPt.x -= c_halfWin_x;
-    prevPt.y -= c_halfWin_y;
+    prevPt -= c_halfWin;
 
     // extract the patch from the first image, compute covariation matrix of derivatives
 
@@ -375,34 +450,68 @@ __kernel void lkSparse_C1_D5(image2d_t I, image2d_t J,
     float A12 = 0;
     float A22 = 0;
 
-    float I_patch[21][21];
-    float dIdx_patch[21][21];
-    float dIdy_patch[21][21];
+    float I_patch[3][3];
+    float dIdx_patch[3][3];
+    float dIdy_patch[3][3];
 
-    for (int yBase = get_local_id(1), i = 0; yBase < c_winSize_y; yBase += get_local_size(1), ++i)
-    {
-        for (int xBase = get_local_id(0), j = 0; xBase < c_winSize_x; xBase += get_local_size(0), ++j)
+        yBase=yid;
         {
-            float x = (prevPt.x + xBase + 0.5f);
-            float y = (prevPt.y + yBase + 0.5f);
-
-            I_patch[i][j] = read_imagef(I, sampler, (float2)(x, y)).x;
-
-            float dIdx = 3.0f * read_imagef(I, sampler, (float2)(x + 1, y - 1)).x + 10.0f * read_imagef(I, sampler, (float2)(x + 1, y)).x + 3.0f * read_imagef(I, sampler, (float2)(x + 1, y + 1)).x -
-                             (3.0f * read_imagef(I, sampler, (float2)(x - 1, y - 1)).x + 10.0f * read_imagef(I, sampler, (float2)(x - 1, y)).x + 3.0f * read_imagef(I, sampler, (float2)(x - 1, y + 1)).x);
-
-            float dIdy = 3.0f * read_imagef(I, sampler, (float2)(x - 1, y + 1)).x + 10.0f * read_imagef(I, sampler, (float2)(x, y + 1)).x + 3.0f * read_imagef(I, sampler, (float2)(x + 1, y + 1)).x -
-                            (3.0f * read_imagef(I, sampler, (float2)(x - 1, y - 1)).x + 10.0f * read_imagef(I, sampler, (float2)(x, y - 1)).x + 3.0f * read_imagef(I, sampler, (float2)(x + 1, y - 1)).x);
-
-            dIdx_patch[i][j] = dIdx;
-            dIdy_patch[i][j] = dIdy;
-
-            A11 += dIdx * dIdx;
-            A12 += dIdx * dIdy;
-            A22 += dIdy * dIdy;
+                xBase=xid;
+                SetPatch(I, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                        &I_patch[0][0], &dIdx_patch[0][0], &dIdy_patch[0][0],
+                                        &A11, &A12, &A22);
+
+
+                xBase+=xsize;
+                SetPatch(I, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                        &I_patch[0][1], &dIdx_patch[0][1], &dIdy_patch[0][1],
+                                        &A11, &A12, &A22);
+
+                xBase+=xsize;
+                if(xBase<c_winSize_x)
+                SetPatch(I, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                        &I_patch[0][2], &dIdx_patch[0][2], &dIdy_patch[0][2],
+                                        &A11, &A12, &A22);
+        }
+        yBase+=ysize;
+        {
+                xBase=xid;
+                SetPatch(I, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                        &I_patch[1][0], &dIdx_patch[1][0], &dIdy_patch[1][0],
+                                        &A11, &A12, &A22);
+
+
+                xBase+=xsize;
+                SetPatch(I, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                        &I_patch[1][1], &dIdx_patch[1][1], &dIdy_patch[1][1],
+                                        &A11, &A12, &A22);
+
+                xBase+=xsize;
+                if(xBase<c_winSize_x)
+                SetPatch(I, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                        &I_patch[1][2], &dIdx_patch[1][2], &dIdy_patch[1][2],
+                                        &A11, &A12, &A22);
+        }
+        yBase+=ysize;
+        if(yBase<c_winSize_y)
+        {
+                xBase=xid;
+                SetPatch(I, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                        &I_patch[2][0], &dIdx_patch[2][0], &dIdy_patch[2][0],
+                                        &A11, &A12, &A22);
+
+
+                xBase+=xsize;
+                SetPatch(I, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                        &I_patch[2][1], &dIdx_patch[2][1], &dIdy_patch[2][1],
+                                        &A11, &A12, &A22);
+
+                xBase+=xsize;
+                if(xBase<c_winSize_x)
+                SetPatch(I, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                        &I_patch[2][2], &dIdx_patch[2][2], &dIdy_patch[2][2],
+                                        &A11, &A12, &A22);
         }
-    }
-
     reduce3(A11, A12, A22, smem1, smem2, smem3, tid);
     barrier(CLK_LOCAL_MEM_FENCE);
 
@@ -412,58 +521,90 @@ __kernel void lkSparse_C1_D5(image2d_t I, image2d_t J,
 
     float D = A11 * A22 - A12 * A12;
 
-    //if (calcErr && GET_MIN_EIGENVALS && tid == 0)
-    //    err[get_group_id(0)] = minEig;
-
     if (D < 1.192092896e-07f)
     {
-        if (level == 0 && tid == 0)
-            status[get_group_id(0)] = 0;
+        if (tid == 0 && level == 0)
+            status[gid] = 0;
 
         return;
     }
 
-    D = 1.f / D;
-
-    A11 *= D;
-    A12 *= D;
-    A22 *= D;
-
-    float2 nextPt = nextPts[get_group_id(0)];
-    nextPt.x *= 2.0f;
-    nextPt.y *= 2.0f;
+    A11 /= D;
+    A12 /= D;
+    A22 /= D;
 
-    nextPt.x -= c_halfWin_x;
-    nextPt.y -= c_halfWin_y;
+    prevPt = nextPts[gid] * 2.0f - c_halfWin;
 
-    for (int k = 0; k < c_iters; ++k)
+    for (k = 0; k < c_iters; ++k)
     {
-        if (nextPt.x < -c_halfWin_x || nextPt.x >= cols || nextPt.y < -c_halfWin_y || nextPt.y >= rows)
+        if (prevPt.x < -c_halfWin.x || prevPt.x >= cols || prevPt.y < -c_halfWin.y || prevPt.y >= rows)
         {
             if (tid == 0 && level == 0)
-                status[get_group_id(0)] = 0;
+                status[gid] = 0;
             return;
         }
 
         float b1 = 0;
         float b2 = 0;
 
-        for (int y = get_local_id(1), i = 0; y < c_winSize_y; y += get_local_size(1), ++i)
-        {
-            for (int x = get_local_id(0), j = 0; x < c_winSize_x; x += get_local_size(0), ++j)
-            {
-                float a = (nextPt.x + x + 0.5f);
-                float b = (nextPt.y + y + 0.5f);
-
-                float I_val = I_patch[i][j];
-                float J_val = read_imagef(J, sampler, (float2)(a, b)).x;
-
-                float diff = (J_val - I_val) * 32.0f;
-
-                b1 += diff * dIdx_patch[i][j];
-                b2 += diff * dIdy_patch[i][j];
-            }
-        }
+                yBase=yid;
+                {
+                        xBase=xid;
+                        GetPatch(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[0][0], &dIdx_patch[0][0], &dIdy_patch[0][0],
+                                                &b1, &b2);
+
+
+                        xBase+=xsize;
+                        GetPatch(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[0][1], &dIdx_patch[0][1], &dIdy_patch[0][1],
+                                                &b1, &b2);
+
+                        xBase+=xsize;
+                        if(xBase<c_winSize_x)
+                        GetPatch(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[0][2], &dIdx_patch[0][2], &dIdy_patch[0][2],
+                                                &b1, &b2);
+                }
+                yBase+=ysize;
+                {
+                        xBase=xid;
+                        GetPatch(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[1][0], &dIdx_patch[1][0], &dIdy_patch[1][0],
+                                                &b1, &b2);
+
+
+                        xBase+=xsize;
+                        GetPatch(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[1][1], &dIdx_patch[1][1], &dIdy_patch[1][1],
+                                                &b1, &b2);
+
+                        xBase+=xsize;
+                        if(xBase<c_winSize_x)
+                        GetPatch(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[1][2], &dIdx_patch[1][2], &dIdy_patch[1][2],
+                                                &b1, &b2);
+                }
+                yBase+=ysize;
+                if(yBase<c_winSize_y)
+                {
+                        xBase=xid;
+                        GetPatch(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[2][0], &dIdx_patch[2][0], &dIdy_patch[2][0],
+                                                &b1, &b2);
+
+
+                        xBase+=xsize;
+                        GetPatch(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[2][1], &dIdx_patch[2][1], &dIdy_patch[2][1],
+                                                &b1, &b2);
+
+                        xBase+=xsize;
+                        if(xBase<c_winSize_x)
+                        GetPatch(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[2][2], &dIdx_patch[2][2], &dIdy_patch[2][2],
+                                                &b1, &b2);
+                }
 
         reduce2(b1, b2, smem1, smem2, tid);
         barrier(CLK_LOCAL_MEM_FENCE);
@@ -475,77 +616,112 @@ __kernel void lkSparse_C1_D5(image2d_t I, image2d_t J,
         delta.x = A12 * b2 - A22 * b1;
         delta.y = A12 * b1 - A11 * b2;
 
-        nextPt.x += delta.x;
-        nextPt.y += delta.y;
+                prevPt += delta;
 
-        if (fabs(delta.x) < 0.01f && fabs(delta.y) < 0.01f)
+        if (fabs(delta.x) < THRESHOLD && fabs(delta.y) < THRESHOLD)
             break;
     }
 
-    float errval = 0.0f;
+    D = 0.0f;
     if (calcErr)
     {
-        for (int y = get_local_id(1), i = 0; y < c_winSize_y; y += get_local_size(1), ++i)
-        {
-            for (int x = get_local_id(0), j = 0; x < c_winSize_x; x += get_local_size(0), ++j)
-            {
-                float a = (nextPt.x + x + 0.5f);
-                float b = (nextPt.y + y + 0.5f);
-
-                float I_val = I_patch[i][j];
-                float J_val = read_imagef(J, sampler, (float2)(a, b)).x;
-
-                float diff = J_val - I_val;
-
-                errval += fabs((float)diff);
-            }
-        }
-
-        reduce1(errval, smem1, tid);
+                yBase=yid;
+                {
+                        xBase=xid;
+                        GetError(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[0][0], &D);
+
+
+                        xBase+=xsize;
+                        GetError(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[0][1], &D);
+
+                        xBase+=xsize;
+                        if(xBase<c_winSize_x)
+                        GetError(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[0][2], &D);
+                }
+                yBase+=ysize;
+                {
+                        xBase=xid;
+                        GetError(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[1][0], &D);
+
+
+                        xBase+=xsize;
+                        GetError(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[1][1], &D);
+
+                        xBase+=xsize;
+                        if(xBase<c_winSize_x)
+                        GetError(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[1][2], &D);
+                }
+                yBase+=ysize;
+                if(yBase<c_winSize_y)
+                {
+                        xBase=xid;
+                        GetError(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[2][0], &D);
+
+
+                        xBase+=xsize;
+                        GetError(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[2][1], &D);
+
+                        xBase+=xsize;
+                        if(xBase<c_winSize_x)
+                        GetError(J, prevPt.x + xBase + 0.5f, prevPt.y + yBase + 0.5f,
+                                                &I_patch[2][2], &D);
+                }
+
+        reduce1(D, smem1, tid);
     }
 
     if (tid == 0)
     {
-        nextPt.x += c_halfWin_x;
-        nextPt.y += c_halfWin_y;
+                prevPt += c_halfWin;
 
-        nextPts[get_group_id(0)] = nextPt;
+        nextPts[gid] = prevPt;
 
-        //if (calcErr && !GET_MIN_EIGENVALS)
-        //    err[get_group_id(0)] = errval;
+        if (calcErr)
+            err[gid] = smem1[0] / (c_winSize_x * c_winSize_y);
     }
+
 }
+
 __kernel void lkSparse_C4_D5(image2d_t I, image2d_t J,
-    __global const float2* prevPts, int prevPtsStep, __global float2* nextPts, int nextPtsStep, __global uchar* status/*, __global float* err*/, const int level, const int rows, const int cols, int PATCH_X, int PATCH_Y, int cn, int c_winSize_x, int c_winSize_y, int c_iters, char calcErr, char GET_MIN_EIGENVALS)
+    __global const float2* prevPts, int prevPtsStep, __global float2* nextPts, int nextPtsStep, __global uchar* status, __global float* err,
+        const int level, const int rows, const int cols, int PATCH_X, int PATCH_Y, int cn, int c_winSize_x, int c_winSize_y, int c_iters, char calcErr)
 {
-    __local float smem1[256];
-    __local float smem2[256];
-    __local float smem3[256];
+    __local float smem1[BUFFER];
+    __local float smem2[BUFFER];
+    __local float smem3[BUFFER];
 
-    int c_halfWin_x = (c_winSize_x - 1) / 2;
-    int c_halfWin_y = (c_winSize_y - 1) / 2;
+        unsigned int xid=get_local_id(0);
+        unsigned int yid=get_local_id(1);
+        unsigned int gid=get_group_id(0);
+        unsigned int xsize=get_local_size(0);
+        unsigned int ysize=get_local_size(1);
+        int xBase, yBase, i, j, k;
 
-    const int tid = get_local_id(1) * get_local_size(0) + get_local_id(0);
+        float2 c_halfWin = (float2)((c_winSize_x - 1)>>1, (c_winSize_y - 1)>>1);
 
-    float2 prevPt = prevPts[get_group_id(0)];
-    prevPt.x *= (1.0f / (1 << level));
-    prevPt.y *= (1.0f / (1 << level));
+    const int tid = mad24(yid, xsize, xid);
 
-    if (prevPt.x < 0 || prevPt.x >= cols || prevPt.y < 0 || prevPt.y >= rows)
+    float2 nextPt = prevPts[gid]/(1<<level);
+
+    if (nextPt.x < 0 || nextPt.x >= cols || nextPt.y < 0 || nextPt.y >= rows)
     {
-        if (level == 0 && tid == 0)
+        if (tid == 0 && level == 0)
         {
-            status[get_group_id(0)] = 0;
-
-            //if (calcErr)
-            //    err[get_group_id(0)] = 0;
+            status[gid] = 0;
         }
 
         return;
     }
 
-    prevPt.x -= c_halfWin_x;
-    prevPt.y -= c_halfWin_y;
+        nextPt -= c_halfWin;
 
     // extract the patch from the first image, compute covariation matrix of derivatives
 
@@ -553,33 +729,70 @@ __kernel void lkSparse_C4_D5(image2d_t I, image2d_t J,
     float A12 = 0;
     float A22 = 0;
 
-    float4 I_patch[21][21];
-    float4 dIdx_patch[21][21];
-    float4 dIdy_patch[21][21];
+    float4 I_patch[8];
+    float4 dIdx_patch[8];
+    float4 dIdy_patch[8];
+        float4 I_add,Dx_add,Dy_add;
 
-    for (int yBase = get_local_id(1), i = 0; yBase < c_winSize_y; yBase += get_local_size(1), ++i)
-    {
-        for (int xBase = get_local_id(0), j = 0; xBase < c_winSize_x; xBase += get_local_size(0), ++j)
+        yBase=yid;
         {
-            float x = (prevPt.x + xBase + 0.5f);
-            float y = (prevPt.y + yBase + 0.5f);
+                xBase=xid;
+                SetPatch4(I, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                        &I_patch[0], &dIdx_patch[0], &dIdy_patch[0],
+                                        &A11, &A12, &A22);
 
-            I_patch[i][j] = read_imagef(I, sampler, (float2)(x, y)).x;
-
-            float4 dIdx = 3.0f * read_imagef(I, sampler, (float2)(x + 1, y - 1)).x + 10.0f * read_imagef(I, sampler, (float2)(x + 1, y)).x + 3.0f * read_imagef(I, sampler, (float2)(x + 1, y + 1)).x -
-                             (3.0f * read_imagef(I, sampler, (float2)(x - 1, y - 1)).x + 10.0f * read_imagef(I, sampler, (float2)(x - 1, y)).x + 3.0f * read_imagef(I, sampler, (float2)(x - 1, y + 1)).x);
 
-            float4 dIdy = 3.0f * read_imagef(I, sampler, (float2)(x - 1, y + 1)).x + 10.0f * read_imagef(I, sampler, (float2)(x, y + 1)).x + 3.0f * read_imagef(I, sampler, (float2)(x + 1, y + 1)).x -
-                            (3.0f * read_imagef(I, sampler, (float2)(x - 1, y - 1)).x + 10.0f * read_imagef(I, sampler, (float2)(x, y - 1)).x + 3.0f * read_imagef(I, sampler, (float2)(x + 1, y - 1)).x);
+                xBase+=xsize;
+                SetPatch4(I, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                        &I_patch[1], &dIdx_patch[1], &dIdy_patch[1],
+                                        &A11, &A12, &A22);
 
-            dIdx_patch[i][j] = dIdx;
-            dIdy_patch[i][j] = dIdy;
+                xBase+=xsize;
+                if(xBase<c_winSize_x)
+                SetPatch4(I, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                        &I_patch[2], &dIdx_patch[2], &dIdy_patch[2],
+                                        &A11, &A12, &A22);
 
-            A11 += (dIdx * dIdx).x + (dIdx * dIdx).y + (dIdx * dIdx).z;
-            A12 += (dIdx * dIdy).x + (dIdx * dIdy).y + (dIdx * dIdy).z;
-            A22 += (dIdy * dIdy).x + (dIdy * dIdy).y + (dIdy * dIdy).z;
         }
-    }
+        yBase+=ysize;
+        {
+                xBase=xid;
+                SetPatch4(I, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                        &I_patch[3], &dIdx_patch[3], &dIdy_patch[3],
+                                        &A11, &A12, &A22);
+
+
+                xBase+=xsize;
+                SetPatch4(I, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                        &I_patch[4], &dIdx_patch[4], &dIdy_patch[4],
+                                        &A11, &A12, &A22);
+
+                xBase+=xsize;
+                if(xBase<c_winSize_x)
+                SetPatch4(I, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                        &I_patch[5], &dIdx_patch[5], &dIdy_patch[5],
+                                        &A11, &A12, &A22);
+        }
+        yBase+=ysize;
+        if(yBase<c_winSize_y)
+        {
+                xBase=xid;
+                SetPatch4(I, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                        &I_patch[6], &dIdx_patch[6], &dIdy_patch[6],
+                                        &A11, &A12, &A22);
+
+
+                xBase+=xsize;
+                SetPatch4(I, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                        &I_patch[7], &dIdx_patch[7], &dIdy_patch[7],
+                                        &A11, &A12, &A22);
+
+                xBase+=xsize;
+                if(xBase<c_winSize_x)
+                SetPatch4(I, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                        &I_add, &Dx_add, &Dy_add,
+                                        &A11, &A12, &A22);
+        }
 
     reduce3(A11, A12, A22, smem1, smem2, smem3, tid);
     barrier(CLK_LOCAL_MEM_FENCE);
@@ -590,58 +803,91 @@ __kernel void lkSparse_C4_D5(image2d_t I, image2d_t J,
 
     float D = A11 * A22 - A12 * A12;
 
-    //if (calcErr && GET_MIN_EIGENVALS && tid == 0)
-    //    err[get_group_id(0)] = minEig;
-
     if (D < 1.192092896e-07f)
     {
-        if (level == 0 && tid == 0)
-            status[get_group_id(0)] = 0;
+        if (tid == 0 && level == 0)
+            status[gid] = 0;
 
         return;
     }
 
-    D = 1.f / D;
-
-    A11 *= D;
-    A12 *= D;
-    A22 *= D;
-
-    float2 nextPt = nextPts[get_group_id(0)];
-    nextPt.x *= 2.0f;
-    nextPt.y *= 2.0f;
+    A11 /= D;
+    A12 /= D;
+    A22 /= D;
 
-    nextPt.x -= c_halfWin_x;
-    nextPt.y -= c_halfWin_y;
+        nextPt = nextPts[gid] * 2.0f - c_halfWin;
 
-    for (int k = 0; k < c_iters; ++k)
+    for (k = 0; k < c_iters; ++k)
     {
-        if (nextPt.x < -c_halfWin_x || nextPt.x >= cols || nextPt.y < -c_halfWin_y || nextPt.y >= rows)
+        if (nextPt.x < -c_halfWin.x || nextPt.x >= cols || nextPt.y < -c_halfWin.y || nextPt.y >= rows)
         {
             if (tid == 0 && level == 0)
-                status[get_group_id(0)] = 0;
+                status[gid] = 0;
             return;
         }
 
         float b1 = 0;
         float b2 = 0;
 
-        for (int y = get_local_id(1), i = 0; y < c_winSize_y; y += get_local_size(1), ++i)
-        {
-            for (int x = get_local_id(0), j = 0; x < c_winSize_x; x += get_local_size(0), ++j)
-            {
-                float a = (nextPt.x + x + 0.5f);
-                float b = (nextPt.y + y + 0.5f);
-
-                float4 I_val = I_patch[i][j];
-                float4 J_val = read_imagef(J, sampler, (float2)(a, b)).x;
-
-                float4 diff = (J_val - I_val) * 32.0f;
+                yBase=yid;
+                {
+                        xBase=xid;
+                        GetPatch4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[0], &dIdx_patch[0], &dIdy_patch[0],
+                                                &b1, &b2);
+
+
+                        xBase+=xsize;
+                        GetPatch4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[1], &dIdx_patch[1], &dIdy_patch[1],
+                                                &b1, &b2);
+
+                        xBase+=xsize;
+                        if(xBase<c_winSize_x)
+                        GetPatch4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[2], &dIdx_patch[2], &dIdy_patch[2],
+                                                &b1, &b2);
+                }
+                yBase+=ysize;
+                {
+                        xBase=xid;
+                        GetPatch4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[3], &dIdx_patch[3], &dIdy_patch[3],
+                                                &b1, &b2);
+
+
+                        xBase+=xsize;
+                        GetPatch4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[4], &dIdx_patch[4], &dIdy_patch[4],
+                                                &b1, &b2);
+
+                        xBase+=xsize;
+                        if(xBase<c_winSize_x)
+                        GetPatch4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[5], &dIdx_patch[5], &dIdy_patch[5],
+                                                &b1, &b2);
+                }
+                yBase+=ysize;
+                if(yBase<c_winSize_y)
+                {
+                        xBase=xid;
+                        GetPatch4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[6], &dIdx_patch[6], &dIdy_patch[6],
+                                                &b1, &b2);
+
+
+                        xBase+=xsize;
+                        GetPatch4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[7], &dIdx_patch[7], &dIdy_patch[7],
+                                                &b1, &b2);
+
+                        xBase+=xsize;
+                        if(xBase<c_winSize_x)
+                        GetPatch4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_add, &Dx_add, &Dy_add,
+                                                &b1, &b2);
+                }
 
-                b1 += (diff * dIdx_patch[i][j]).x + (diff * dIdx_patch[i][j]).y + (diff * dIdx_patch[i][j]).z;
-                b2 += (diff * dIdy_patch[i][j]).x + (diff * dIdy_patch[i][j]).y + (diff * dIdy_patch[i][j]).z;
-            }
-        }
 
         reduce2(b1, b2, smem1, smem2, tid);
         barrier(CLK_LOCAL_MEM_FENCE);
@@ -653,52 +899,83 @@ __kernel void lkSparse_C4_D5(image2d_t I, image2d_t J,
         delta.x = A12 * b2 - A22 * b1;
         delta.y = A12 * b1 - A11 * b2;
 
-        nextPt.x += delta.x;
-        nextPt.y += delta.y;
+                nextPt +=delta;
 
-        if (fabs(delta.x) < 0.01f && fabs(delta.y) < 0.01f)
+        if (fabs(delta.x) < THRESHOLD && fabs(delta.y) < THRESHOLD)
             break;
     }
 
-    float errval = 0.0f;
+    D = 0.0f;
     if (calcErr)
     {
-        for (int y = get_local_id(1), i = 0; y < c_winSize_y; y += get_local_size(1), ++i)
-        {
-            for (int x = get_local_id(0), j = 0; x < c_winSize_x; x += get_local_size(0), ++j)
-            {
-                float a = (nextPt.x + x + 0.5f);
-                float b = (nextPt.y + y + 0.5f);
-
-                float4 I_val = I_patch[i][j];
-                float4 J_val = read_imagef(J, sampler, (float2)(a, b)).x;
-
-                float4 diff = J_val - I_val;
-
-                errval += fabs(diff.x) + fabs(diff.y) + fabs(diff.z);
-            }
-        }
-
-        reduce1(errval, smem1, tid);
+                yBase=yid;
+                {
+                        xBase=xid;
+                        GetError4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[0], &D);
+
+
+                        xBase+=xsize;
+                        GetError4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[1], &D);
+
+                        xBase+=xsize;
+                        if(xBase<c_winSize_x)
+                        GetError4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[2], &D);
+                }
+                yBase+=ysize;
+                {
+                        xBase=xid;
+                        GetError4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[3], &D);
+
+
+                        xBase+=xsize;
+                        GetError4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[4], &D);
+
+                        xBase+=xsize;
+                        if(xBase<c_winSize_x)
+                        GetError4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[5], &D);
+                }
+                yBase+=ysize;
+                if(yBase<c_winSize_y)
+                {
+                        xBase=xid;
+                        GetError4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[6], &D);
+
+
+                        xBase+=xsize;
+                        GetError4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_patch[7], &D);
+
+                        xBase+=xsize;
+                        if(xBase<c_winSize_x)
+                        GetError4(J, nextPt.x + xBase + 0.5f, nextPt.y + yBase + 0.5f,
+                                                &I_add, &D);
+                }
+
+        reduce1(D, smem1, tid);
     }
 
     if (tid == 0)
     {
-        nextPt.x += c_halfWin_x;
-        nextPt.y += c_halfWin_y;
-
-        nextPts[get_group_id(0)] = nextPt;
+                nextPt += c_halfWin;
+        nextPts[gid] = nextPt;
 
-        //if (calcErr && !GET_MIN_EIGENVALS)
-        //    err[get_group_id(0)] = errval;
+        if (calcErr)
+            err[gid] = smem1[0] / (3 * c_winSize_x * c_winSize_y);
     }
 }
 
 __kernel void lkDense_C1_D0(image2d_t I, image2d_t J, __global float* u, int uStep, __global float* v, int vStep, __global const float* prevU, int prevUStep, __global const float* prevV, int prevVStep,
     const int rows, const int cols, /*__global float* err, int errStep, int cn,*/ int c_winSize_x, int c_winSize_y, int c_iters, char calcErr)
 {
-    int c_halfWin_x = (c_winSize_x - 1) / 2;
-    int c_halfWin_y = (c_winSize_y - 1) / 2;
+        int c_halfWin_x = (c_winSize_x - 1) / 2;
+        int c_halfWin_y = (c_winSize_y - 1) / 2;
 
     const int patchWidth  = get_local_size(0) + 2 * c_halfWin_x;
     const int patchHeight = get_local_size(1) + 2 * c_halfWin_y;
@@ -712,7 +989,7 @@ __kernel void lkDense_C1_D0(image2d_t I, image2d_t J, __global float* u, int uSt
     const int xBase = get_group_id(0) * get_local_size(0);
     const int yBase = get_group_id(1) * get_local_size(1);
 
-    sampler_t sampleri    = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST;
+        sampler_t sampleri    = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST;
 
     for (int i = get_local_id(1); i < patchHeight; i += get_local_size(1))
     {
diff --git a/modules/ocl/src/kernels/pyrlk_no_image.cl b/modules/ocl/src/kernels/pyrlk_no_image.cl
new file mode 100644 (file)
index 0000000..98a11b5
--- /dev/null
@@ -0,0 +1,764 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2010-2012, Multicoreware, Inc., all rights reserved.
+// Copyright (C) 2010-2012, Advanced Micro Devices, Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// @Authors
+//    Sen Liu, sen@multicorewareinc.com
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other oclMaterials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors as is and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#define        BUFFER  256
+void reduce3(float val1, float val2, float val3, __local float *smem1, __local float *smem2, __local float *smem3, int tid)
+{
+        smem1[tid] = val1;
+        smem2[tid] = val2;
+        smem3[tid] = val3;
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+#if    BUFFER > 128
+
+        if (tid < 128)
+        {
+                smem1[tid] = val1 += smem1[tid + 128];
+                smem2[tid] = val2 += smem2[tid + 128];
+                smem3[tid] = val3 += smem3[tid + 128];
+        }
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+#endif
+
+#if    BUFFER > 64
+
+        if (tid < 64)
+        {
+                smem1[tid] = val1 += smem1[tid + 64];
+                smem2[tid] = val2 += smem2[tid + 64];
+                smem3[tid] = val3 += smem3[tid + 64];
+        }
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+#endif
+
+        if (tid < 32)
+        {
+                smem1[tid] = val1 += smem1[tid + 32];
+                smem2[tid] = val2 += smem2[tid + 32];
+                smem3[tid] = val3 += smem3[tid + 32];
+        }
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+        if (tid < 16)
+        {
+                smem1[tid] = val1 += smem1[tid + 16];
+                smem2[tid] = val2 += smem2[tid + 16];
+                smem3[tid] = val3 += smem3[tid + 16];
+        }
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+        if (tid < 8)
+        {
+                volatile __local float *vmem1 = smem1;
+                volatile __local float *vmem2 = smem2;
+                volatile __local float *vmem3 = smem3;
+
+                vmem1[tid] = val1 += vmem1[tid + 8];
+                vmem2[tid] = val2 += vmem2[tid + 8];
+                vmem3[tid] = val3 += vmem3[tid + 8];
+
+                vmem1[tid] = val1 += vmem1[tid + 4];
+                vmem2[tid] = val2 += vmem2[tid + 4];
+                vmem3[tid] = val3 += vmem3[tid + 4];
+
+                vmem1[tid] = val1 += vmem1[tid + 2];
+                vmem2[tid] = val2 += vmem2[tid + 2];
+                vmem3[tid] = val3 += vmem3[tid + 2];
+
+                vmem1[tid] = val1 += vmem1[tid + 1];
+                vmem2[tid] = val2 += vmem2[tid + 1];
+                vmem3[tid] = val3 += vmem3[tid + 1];
+        }
+}
+
+void reduce2(float val1, float val2, __local float *smem1, __local float *smem2, int tid)
+{
+        smem1[tid] = val1;
+        smem2[tid] = val2;
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+#if    BUFFER > 128
+
+        if (tid < 128)
+        {
+                smem1[tid] = val1 += smem1[tid + 128];
+                smem2[tid] = val2 += smem2[tid + 128];
+        }
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+#endif
+
+#if    BUFFER > 64
+
+        if (tid < 64)
+        {
+                smem1[tid] = val1 += smem1[tid + 64];
+                smem2[tid] = val2 += smem2[tid + 64];
+        }
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+#endif
+
+        if (tid < 32)
+        {
+                smem1[tid] = val1 += smem1[tid + 32];
+                smem2[tid] = val2 += smem2[tid + 32];
+        }
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+        if (tid < 16)
+        {
+                smem1[tid] = val1 += smem1[tid + 16];
+                smem2[tid] = val2 += smem2[tid + 16];
+        }
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+        if (tid < 8)
+        {
+                volatile __local float *vmem1 = smem1;
+                volatile __local float *vmem2 = smem2;
+
+                vmem1[tid] = val1 += vmem1[tid + 8];
+                vmem2[tid] = val2 += vmem2[tid + 8];
+
+                vmem1[tid] = val1 += vmem1[tid + 4];
+                vmem2[tid] = val2 += vmem2[tid + 4];
+
+                vmem1[tid] = val1 += vmem1[tid + 2];
+                vmem2[tid] = val2 += vmem2[tid + 2];
+
+                vmem1[tid] = val1 += vmem1[tid + 1];
+                vmem2[tid] = val2 += vmem2[tid + 1];
+        }
+}
+
+void reduce1(float val1, __local float *smem1, int tid)
+{
+        smem1[tid] = val1;
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+#if    BUFFER > 128
+
+        if (tid < 128)
+        {
+                smem1[tid] = val1 += smem1[tid + 128];
+        }
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+#endif
+
+#if    BUFFER > 64
+
+        if (tid < 64)
+        {
+                smem1[tid] = val1 += smem1[tid + 64];
+        }
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+#endif
+
+        if (tid < 32)
+        {
+                smem1[tid] = val1 += smem1[tid + 32];
+        }
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+        if (tid < 16)
+        {
+                volatile __local float *vmem1 = smem1;
+
+                vmem1[tid] = val1 += vmem1[tid + 16];
+        }
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+        if (tid < 8)
+        {
+                volatile __local float *vmem1 = smem1;
+
+                vmem1[tid] = val1 += vmem1[tid + 8];
+                vmem1[tid] = val1 += vmem1[tid + 4];
+                vmem1[tid] = val1 += vmem1[tid + 2];
+                vmem1[tid] = val1 += vmem1[tid + 1];
+        }
+}
+
+#define SCALE (1.0f / (1 << 20))
+#define        THRESHOLD       0.01f
+#define        DIMENSION       21
+
+float readImage2Df_C1(__global const float *image,  const float x,  const float y,  const int rows,  const int cols, const int elemCntPerRow)
+{
+        float2 coor = (float2)(x, y);
+
+        int i0 = clamp((int)floor(coor.x), 0, cols - 1);
+        int j0 = clamp((int)floor(coor.y), 0, rows - 1);
+        int i1 = clamp((int)floor(coor.x) + 1, 0, cols - 1);
+        int j1 = clamp((int)floor(coor.y) + 1, 0, rows - 1);
+        float a = coor.x - floor(coor.x);
+        float b = coor.y - floor(coor.y);
+
+        return (1 - a) * (1 - b) * image[mad24(j0, elemCntPerRow, i0)]
+               + a * (1 - b) * image[mad24(j0, elemCntPerRow, i1)]
+               + (1 - a) * b * image[mad24(j1, elemCntPerRow, i0)]
+               + a * b * image[mad24(j1, elemCntPerRow, i1)];
+}
+
+__kernel void lkSparse_C1_D5(__global const float *I, __global const float *J,
+                             __global const float2 *prevPts, int prevPtsStep, __global float2 *nextPts, int nextPtsStep, __global uchar *status, __global float *err,
+                             const int level, const int rows, const int cols, const int elemCntPerRow,
+                             int PATCH_X, int PATCH_Y, int cn, int c_winSize_x, int c_winSize_y, int c_iters, char calcErr)
+{
+        __local float smem1[BUFFER];
+        __local float smem2[BUFFER];
+        __local float smem3[BUFFER];
+
+        float2 c_halfWin = (float2)((c_winSize_x - 1) >> 1, (c_winSize_y - 1) >> 1);
+
+        const int tid = mad24(get_local_id(1), get_local_size(0), get_local_id(0));
+
+        float2 prevPt = prevPts[get_group_id(0)] * (1.0f / (1 << level));
+
+        if (prevPt.x < 0 || prevPt.x >= cols || prevPt.y < 0 || prevPt.y >= rows)
+        {
+                if (tid == 0 && level == 0)
+                {
+                        status[get_group_id(0)] = 0;
+                }
+
+                return;
+        }
+
+        prevPt -= c_halfWin;
+
+        // extract the patch from the first image, compute covariation matrix of derivatives
+
+        float A11 = 0;
+        float A12 = 0;
+        float A22 = 0;
+
+        float I_patch[1][3];
+        float dIdx_patch[1][3];
+        float dIdy_patch[1][3];
+
+        for (int yBase = get_local_id(1), i = 0; yBase < c_winSize_y; yBase += get_local_size(1), ++i)
+        {
+                for (int xBase = get_local_id(0), j = 0; xBase < c_winSize_x; xBase += get_local_size(0), ++j)
+                {
+                        float x = (prevPt.x + xBase);
+                        float y = (prevPt.y + yBase);
+
+                        I_patch[i][j] = readImage2Df_C1(I, x, y, rows, cols, elemCntPerRow);
+                        float dIdx = 3.0f * readImage2Df_C1(I, x + 1, y - 1, rows, cols, elemCntPerRow) + 10.0f * readImage2Df_C1(I, x + 1, y, rows, cols, elemCntPerRow) + 3.0f * readImage2Df_C1(I, x + 1, y + 1, rows, cols, elemCntPerRow) -
+                                     (3.0f * readImage2Df_C1(I, x - 1, y - 1, rows, cols, elemCntPerRow) + 10.0f * readImage2Df_C1(I, x - 1, y, rows, cols, elemCntPerRow) + 3.0f * readImage2Df_C1(I, x - 1, y + 1, rows, cols, elemCntPerRow));
+
+                        float dIdy = 3.0f * readImage2Df_C1(I, x - 1, y + 1, rows, cols, elemCntPerRow) + 10.0f * readImage2Df_C1(I, x, y + 1, rows, cols, elemCntPerRow) + 3.0f * readImage2Df_C1(I, x + 1, y + 1, rows, cols, elemCntPerRow) -
+                                     (3.0f * readImage2Df_C1(I, x - 1, y - 1, rows, cols, elemCntPerRow) + 10.0f * readImage2Df_C1(I, x, y - 1, rows, cols, elemCntPerRow) + 3.0f * readImage2Df_C1(I, x + 1, y - 1, rows, cols, elemCntPerRow));
+
+                        dIdx_patch[i][j] = dIdx;
+                        dIdy_patch[i][j] = dIdy;
+
+                        A11 += dIdx * dIdx;
+                        A12 += dIdx * dIdy;
+                        A22 += dIdy * dIdy;
+                }
+        }
+
+        reduce3(A11, A12, A22, smem1, smem2, smem3, tid);
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+        A11 = smem1[0];
+        A12 = smem2[0];
+        A22 = smem3[0];
+
+        float D = A11 * A22 - A12 * A12;
+
+        if (D < 1.192092896e-07f)
+        {
+                if (tid == 0 && level == 0)
+                {
+                        status[get_group_id(0)] = 0;
+                }
+
+                return;
+        }
+
+        D = 1.f / D;
+
+        A11 *= D;
+        A12 *= D;
+        A22 *= D;
+
+        float2 nextPt = nextPts[get_group_id(0)];
+        nextPt = nextPt * 2.0f - c_halfWin;
+
+        for (int k = 0; k < c_iters; ++k)
+        {
+                if (nextPt.x < -c_halfWin.x || nextPt.x >= cols || nextPt.y < -c_halfWin.y || nextPt.y >= rows)
+                {
+                        if (tid == 0 && level == 0)
+                        {
+                                status[get_group_id(0)] = 0;
+                        }
+
+                        return;
+                }
+
+                float b1 = 0;
+                float b2 = 0;
+
+                for (int y = get_local_id(1), i = 0; y < c_winSize_y; y += get_local_size(1), ++i)
+                {
+                        for (int x = get_local_id(0), j = 0; x < c_winSize_x; x += get_local_size(0), ++j)
+                        {
+                                float diff = (readImage2Df_C1(J, nextPt.x + x, nextPt.y + y, rows, cols, elemCntPerRow) - I_patch[i][j]) * 32.0f;
+
+                                b1 += diff * dIdx_patch[i][j];
+                                b2 += diff * dIdy_patch[i][j];
+                        }
+                }
+
+                reduce2(b1, b2, smem1, smem2, tid);
+                barrier(CLK_LOCAL_MEM_FENCE);
+
+                b1 = smem1[0];
+                b2 = smem2[0];
+
+                float2 delta;
+                delta.x = A12 * b2 - A22 * b1;
+                delta.y = A12 * b1 - A11 * b2;
+
+                nextPt += delta;
+
+                //if (fabs(delta.x) < THRESHOLD && fabs(delta.y) < THRESHOLD)
+                //    break;
+        }
+
+        float errval = 0.0f;
+
+        if (calcErr)
+        {
+                for (int y = get_local_id(1), i = 0; y < c_winSize_y; y += get_local_size(1), ++i)
+                {
+                        for (int x = get_local_id(0), j = 0; x < c_winSize_x; x += get_local_size(0), ++j)
+                        {
+                                float diff = readImage2Df_C1(J, nextPt.x + x, nextPt.y + y, rows, cols, elemCntPerRow) - I_patch[i][j];
+
+                                errval += fabs(diff);
+                        }
+                }
+
+                reduce1(errval, smem1, tid);
+        }
+
+        if (tid == 0)
+        {
+                nextPt += c_halfWin;
+
+                nextPts[get_group_id(0)] = nextPt;
+
+                if (calcErr)
+                {
+                        err[get_group_id(0)] = smem1[0] / (c_winSize_x * c_winSize_y);
+                }
+        }
+}
+
+float4 readImage2Df_C4(__global const float4 *image,  const float x,  const float y,  const int rows,  const int cols, const int elemCntPerRow)
+{
+        float2 coor = (float2)(x, y);
+
+        int i0 = clamp((int)floor(coor.x), 0, cols - 1);
+        int j0 = clamp((int)floor(coor.y), 0, rows - 1);
+        int i1 = clamp((int)floor(coor.x) + 1, 0, cols - 1);
+        int j1 = clamp((int)floor(coor.y) + 1, 0, rows - 1);
+        float a = coor.x - floor(coor.x);
+        float b = coor.y - floor(coor.y);
+
+        return (1 - a) * (1 - b) * image[mad24(j0, elemCntPerRow, i0)]
+               + a * (1 - b) * image[mad24(j0, elemCntPerRow, i1)]
+               + (1 - a) * b * image[mad24(j1, elemCntPerRow, i0)]
+               + a * b * image[mad24(j1, elemCntPerRow, i1)];
+}
+
+__kernel void lkSparse_C4_D5(__global const float *I, __global const float *J,
+                             __global const float2 *prevPts, int prevPtsStep, __global float2 *nextPts, int nextPtsStep, __global uchar *status, __global float *err,
+                             const int level, const int rows, const int cols, const int elemCntPerRow,
+                             int PATCH_X, int PATCH_Y, int cn, int c_winSize_x, int c_winSize_y, int c_iters, char calcErr)
+{
+        __local float smem1[BUFFER];
+        __local float smem2[BUFFER];
+        __local float smem3[BUFFER];
+
+        float2 c_halfWin = (float2)((c_winSize_x - 1) >> 1, (c_winSize_y - 1) >> 1);
+
+        const int tid = mad24(get_local_id(1), get_local_size(0), get_local_id(0));
+
+        float2 prevPt = prevPts[get_group_id(0)] * (1.0f / (1 << level));
+
+        if (prevPt.x < 0 || prevPt.x >= cols || prevPt.y < 0 || prevPt.y >= rows)
+        {
+                if (tid == 0 && level == 0)
+                {
+                        status[get_group_id(0)] = 0;
+                }
+
+                return;
+        }
+
+        prevPt -= c_halfWin;
+
+        // extract the patch from the first image, compute covariation matrix of derivatives
+
+        float A11 = 0;
+        float A12 = 0;
+        float A22 = 0;
+
+        float4 I_patch[1][3];
+        float4 dIdx_patch[1][3];
+        float4 dIdy_patch[1][3];
+
+        __global float4 *ptrI = (__global float4 *)I;
+
+        for (int yBase = get_local_id(1), i = 0; yBase < c_winSize_y; yBase += get_local_size(1), ++i)
+        {
+                for (int xBase = get_local_id(0), j = 0; xBase < c_winSize_x; xBase += get_local_size(0), ++j)
+                {
+                        float x = (prevPt.x + xBase);
+                        float y = (prevPt.y + yBase);
+
+                        I_patch[i][j] = readImage2Df_C4(ptrI, x, y, rows, cols, elemCntPerRow);
+
+                        float4 dIdx = 3.0f * readImage2Df_C4(ptrI, x + 1, y - 1, rows, cols, elemCntPerRow) + 10.0f * readImage2Df_C4(ptrI, x + 1, y, rows, cols, elemCntPerRow) + 3.0f * readImage2Df_C4(ptrI, x + 1, y + 1, rows, cols, elemCntPerRow) -
+                                      (3.0f * readImage2Df_C4(ptrI, x - 1, y - 1, rows, cols, elemCntPerRow) + 10.0f * readImage2Df_C4(ptrI, x - 1, y, rows, cols, elemCntPerRow) + 3.0f * readImage2Df_C4(ptrI, x - 1, y + 1, rows, cols, elemCntPerRow));
+
+                        float4 dIdy = 3.0f * readImage2Df_C4(ptrI, x - 1, y + 1, rows, cols, elemCntPerRow) + 10.0f * readImage2Df_C4(ptrI, x, y + 1, rows, cols, elemCntPerRow) + 3.0f * readImage2Df_C4(ptrI, x + 1, y + 1, rows, cols, elemCntPerRow) -
+                                      (3.0f * readImage2Df_C4(ptrI, x - 1, y - 1, rows, cols, elemCntPerRow) + 10.0f * readImage2Df_C4(ptrI, x, y - 1, rows, cols, elemCntPerRow) + 3.0f * readImage2Df_C4(ptrI, x + 1, y - 1, rows, cols, elemCntPerRow));
+
+                        dIdx_patch[i][j] = dIdx;
+                        dIdy_patch[i][j] = dIdy;
+
+                        A11 += (dIdx * dIdx).x + (dIdx * dIdx).y + (dIdx * dIdx).z;
+                        A12 += (dIdx * dIdy).x + (dIdx * dIdy).y + (dIdx * dIdy).z;
+                        A22 += (dIdy * dIdy).x + (dIdy * dIdy).y + (dIdy * dIdy).z;
+                }
+        }
+
+        reduce3(A11, A12, A22, smem1, smem2, smem3, tid);
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+        A11 = smem1[0];
+        A12 = smem2[0];
+        A22 = smem3[0];
+
+        float D = A11 * A22 - A12 * A12;
+        //pD[get_group_id(0)] = D;
+
+        if (D < 1.192092896e-07f)
+        {
+                if (tid == 0 && level == 0)
+                {
+                        status[get_group_id(0)] = 0;
+                }
+
+                return;
+        }
+
+        D = 1.f / D;
+
+        A11 *= D;
+        A12 *= D;
+        A22 *= D;
+
+        float2 nextPt = nextPts[get_group_id(0)];
+
+        nextPt = nextPt * 2.0f - c_halfWin;
+
+        __global float4 *ptrJ = (__global float4 *)J;
+
+        for (int k = 0; k < c_iters; ++k)
+        {
+                if (nextPt.x < -c_halfWin.x || nextPt.x >= cols || nextPt.y < -c_halfWin.y || nextPt.y >= rows)
+                {
+                        if (tid == 0 && level == 0)
+                        {
+                                status[get_group_id(0)] = 0;
+                        }
+
+                        return;
+                }
+
+                float b1 = 0;
+                float b2 = 0;
+
+                for (int y = get_local_id(1), i = 0; y < c_winSize_y; y += get_local_size(1), ++i)
+                {
+                        for (int x = get_local_id(0), j = 0; x < c_winSize_x; x += get_local_size(0), ++j)
+                        {
+                                float4 diff = (readImage2Df_C4(ptrJ, nextPt.x + x, nextPt.y + y, rows, cols, elemCntPerRow) - I_patch[i][j]) * 32.0f;
+
+                                b1 += (diff * dIdx_patch[i][j]).x + (diff * dIdx_patch[i][j]).y + (diff * dIdx_patch[i][j]).z;
+                                b2 += (diff * dIdy_patch[i][j]).x + (diff * dIdy_patch[i][j]).y + (diff * dIdy_patch[i][j]).z;
+                        }
+                }
+
+                reduce2(b1, b2, smem1, smem2, tid);
+                barrier(CLK_LOCAL_MEM_FENCE);
+
+                b1 = smem1[0];
+                b2 = smem2[0];
+
+                float2 delta;
+                delta.x = A12 * b2 - A22 * b1;
+                delta.y = A12 * b1 - A11 * b2;
+
+                nextPt += delta;
+
+                //if (fabs(delta.x) < THRESHOLD && fabs(delta.y) < THRESHOLD)
+                //    break;
+        }
+
+        float errval = 0.0f;
+
+        if (calcErr)
+        {
+                for (int y = get_local_id(1), i = 0; y < c_winSize_y; y += get_local_size(1), ++i)
+                {
+                        for (int x = get_local_id(0), j = 0; x < c_winSize_x; x += get_local_size(0), ++j)
+                        {
+                                float4 diff = readImage2Df_C4(ptrJ, nextPt.x + x, nextPt.y + y, rows, cols, elemCntPerRow) - I_patch[i][j];
+
+                                errval += fabs(diff.x) + fabs(diff.y) + fabs(diff.z);
+                        }
+                }
+
+                reduce1(errval, smem1, tid);
+        }
+
+        if (tid == 0)
+        {
+                nextPt += c_halfWin;
+                nextPts[get_group_id(0)] = nextPt;
+
+                if (calcErr)
+                {
+                        err[get_group_id(0)] = smem1[0] / (3 * c_winSize_x * c_winSize_y);
+                }
+        }
+}
+
+int readImage2Di_C1(__global const int *image, float2 coor,  int2 size, const int elemCntPerRow)
+{
+        int i = clamp((int)floor(coor.x), 0, size.x - 1);
+        int j = clamp((int)floor(coor.y), 0, size.y - 1);
+        return image[mad24(j, elemCntPerRow, i)];
+}
+
+__kernel void lkDense_C1_D0(__global const int *I, __global const int *J, __global float *u, int uStep, __global float *v, int vStep, __global const float *prevU, int prevUStep, __global const float *prevV, int prevVStep,
+                            const int rows, const int cols, /*__global float* err, int errStep, int cn,*/
+                            const int elemCntPerRow, int c_winSize_x, int c_winSize_y, int c_iters, char calcErr)
+{
+        int c_halfWin_x = (c_winSize_x - 1) / 2;
+        int c_halfWin_y = (c_winSize_y - 1) / 2;
+
+        const int patchWidth  = get_local_size(0) + 2 * c_halfWin_x;
+        const int patchHeight = get_local_size(1) + 2 * c_halfWin_y;
+
+        __local int smem[8192];
+
+        __local int *I_patch = smem;
+        __local int *dIdx_patch = I_patch + patchWidth * patchHeight;
+        __local int *dIdy_patch = dIdx_patch + patchWidth * patchHeight;
+
+        const int xBase = get_group_id(0) * get_local_size(0);
+        const int yBase = get_group_id(1) * get_local_size(1);
+        int2 size = (int2)(cols, rows);
+
+        for (int i = get_local_id(1); i < patchHeight; i += get_local_size(1))
+        {
+                for (int j = get_local_id(0); j < patchWidth; j += get_local_size(0))
+                {
+                        float x = xBase - c_halfWin_x + j + 0.5f;
+                        float y = yBase - c_halfWin_y + i + 0.5f;
+
+                        I_patch[i * patchWidth + j] = readImage2Di_C1(I, (float2)(x, y), size, elemCntPerRow);
+
+                        // Sharr Deriv
+
+                        dIdx_patch[i * patchWidth + j] = 3 * readImage2Di_C1(I, (float2)(x + 1, y - 1), size, elemCntPerRow) + 10 * readImage2Di_C1(I, (float2)(x + 1, y), size, elemCntPerRow) + 3 * readImage2Di_C1(I, (float2)(x + 1, y + 1), size, elemCntPerRow) -
+                                                         (3 * readImage2Di_C1(I, (float2)(x - 1, y - 1), size, elemCntPerRow) + 10 * readImage2Di_C1(I, (float2)(x - 1, y), size, elemCntPerRow) + 3 * readImage2Di_C1(I, (float2)(x - 1, y + 1), size, elemCntPerRow));
+
+                        dIdy_patch[i * patchWidth + j] = 3 * readImage2Di_C1(I, (float2)(x - 1, y + 1), size, elemCntPerRow) + 10 * readImage2Di_C1(I, (float2)(x, y + 1), size, elemCntPerRow) + 3 * readImage2Di_C1(I, (float2)(x + 1, y + 1), size, elemCntPerRow) -
+                                                         (3 * readImage2Di_C1(I, (float2)(x - 1, y - 1), size, elemCntPerRow) + 10 * readImage2Di_C1(I, (float2)(x, y - 1), size, elemCntPerRow) + 3 * readImage2Di_C1(I, (float2)(x + 1, y - 1), size, elemCntPerRow));
+                }
+        }
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+        // extract the patch from the first image, compute covariation matrix of derivatives
+
+        const int x = get_global_id(0);
+        const int y = get_global_id(1);
+
+        if (x >= cols || y >= rows)
+        {
+                return;
+        }
+
+        int A11i = 0;
+        int A12i = 0;
+        int A22i = 0;
+
+        for (int i = 0; i < c_winSize_y; ++i)
+        {
+                for (int j = 0; j < c_winSize_x; ++j)
+                {
+                        int dIdx = dIdx_patch[(get_local_id(1) + i) * patchWidth + (get_local_id(0) + j)];
+                        int dIdy = dIdy_patch[(get_local_id(1) + i) * patchWidth + (get_local_id(0) + j)];
+
+                        A11i += dIdx * dIdx;
+                        A12i += dIdx * dIdy;
+                        A22i += dIdy * dIdy;
+                }
+        }
+
+        float A11 = A11i;
+        float A12 = A12i;
+        float A22 = A22i;
+
+        float D = A11 * A22 - A12 * A12;
+
+        //if (calcErr && GET_MIN_EIGENVALS)
+        //    (err + y * errStep)[x] = minEig;
+
+        if (D < 1.192092896e-07f)
+        {
+                //if (calcErr)
+                //    err(y, x) = 3.402823466e+38f;
+
+                return;
+        }
+
+        D = 1.f / D;
+
+        A11 *= D;
+        A12 *= D;
+        A22 *= D;
+
+        float2 nextPt;
+        nextPt.x = x + prevU[y / 2 * prevUStep / 4 + x / 2] * 2.0f;
+        nextPt.y = y + prevV[y / 2 * prevVStep / 4 + x / 2] * 2.0f;
+
+        for (int k = 0; k < c_iters; ++k)
+        {
+                if (nextPt.x < 0 || nextPt.x >= cols || nextPt.y < 0 || nextPt.y >= rows)
+                {
+                        //if (calcErr)
+                        //    err(y, x) = 3.402823466e+38f;
+
+                        return;
+                }
+
+                int b1 = 0;
+                int b2 = 0;
+
+                for (int i = 0; i < c_winSize_y; ++i)
+                {
+                        for (int j = 0; j < c_winSize_x; ++j)
+                        {
+                                int iI = I_patch[(get_local_id(1) + i) * patchWidth + get_local_id(0) + j];
+                                int iJ = readImage2Di_C1(J, (float2)(nextPt.x - c_halfWin_x + j + 0.5f, nextPt.y - c_halfWin_y + i + 0.5f), size, elemCntPerRow);
+
+                                int diff = (iJ - iI) * 32;
+
+                                int dIdx = dIdx_patch[(get_local_id(1) + i) * patchWidth + (get_local_id(0) + j)];
+                                int dIdy = dIdy_patch[(get_local_id(1) + i) * patchWidth + (get_local_id(0) + j)];
+
+                                b1 += diff * dIdx;
+                                b2 += diff * dIdy;
+                        }
+                }
+
+                float2 delta;
+                delta.x = A12 * b2 - A22 * b1;
+                delta.y = A12 * b1 - A11 * b2;
+
+                nextPt.x += delta.x;
+                nextPt.y += delta.y;
+
+                if (fabs(delta.x) < 0.01f && fabs(delta.y) < 0.01f)
+                {
+                        break;
+                }
+        }
+
+        u[y * uStep / 4 + x] = nextPt.x - x;
+        v[y * vStep / 4 + x] = nextPt.y - y;
+
+        if (calcErr)
+        {
+                int errval = 0;
+
+                for (int i = 0; i < c_winSize_y; ++i)
+                {
+                        for (int j = 0; j < c_winSize_x; ++j)
+                        {
+                                int iI = I_patch[(get_local_id(1) + i) * patchWidth + get_local_id(0) + j];
+                                int iJ = readImage2Di_C1(J, (float2)(nextPt.x - c_halfWin_x + j + 0.5f, nextPt.y - c_halfWin_y + i + 0.5f), size, elemCntPerRow);
+
+                                errval += abs(iJ - iI);
+                        }
+                }
+
+                //err[y * errStep / 4 + x] = static_cast<float>(errval) / (c_winSize_x * c_winSize_y);
+        }
+}
index dac303c..7165a8c 100644 (file)
@@ -48,23 +48,24 @@ using namespace cv::ocl;
 
 #if !defined (HAVE_OPENCL)
 
-void cv::ocl::PyrLKOpticalFlow::sparse(const oclMat &, const oclMat &, const oclMat &, oclMat &, oclMat &, oclMat *) {  }
+void cv::ocl::PyrLKOpticalFlow::sparse(const oclMat &, const oclMat &, const oclMat &, oclMat &, oclMat &, oclMat &) {  }
 void cv::ocl::PyrLKOpticalFlow::dense(const oclMat &, const oclMat &, oclMat &, oclMat &, oclMat *) {  }
 
 #else /* !defined (HAVE_OPENCL) */
 
 namespace cv
 {
-    namespace ocl
-    {
-        ///////////////////////////OpenCL kernel strings///////////////////////////
-        extern const char *pyrlk;
-        extern const char *operator_setTo;
-        extern const char *operator_convertTo;
-        extern const char *operator_copyToM;
-        extern const char *arithm_mul;
-        extern const char *pyr_down;
-    }
+namespace ocl
+{
+///////////////////////////OpenCL kernel strings///////////////////////////
+extern const char *pyrlk;
+extern const char *pyrlk_no_image;
+extern const char *operator_setTo;
+extern const char *operator_convertTo;
+extern const char *operator_copyToM;
+extern const char *arithm_mul;
+extern const char *pyr_down;
+}
 }
 
 struct dim3
@@ -84,26 +85,26 @@ struct int2
 
 namespace
 {
-    void calcPatchSize(cv::Size winSize, int cn, dim3 &block, dim3 &patch, bool isDeviceArch11)
-    {
-        winSize.width *= cn;
+void calcPatchSize(cv::Size winSize, int cn, dim3 &block, dim3 &patch, bool isDeviceArch11)
+{
+    winSize.width *= cn;
 
-        if (winSize.width > 32 && winSize.width > 2 * winSize.height)
-        {
-            block.x = isDeviceArch11 ? 16 : 32;
-            block.y = 8;
-        }
-        else
-        {
-            block.x = 16;
-            block.y = isDeviceArch11 ? 8 : 16;
-        }
+    if (winSize.width > 32 && winSize.width > 2 * winSize.height)
+    {
+        block.x = isDeviceArch11 ? 16 : 32;
+        block.y = 8;
+    }
+    else
+    {
+        block.x = 16;
+        block.y = isDeviceArch11 ? 8 : 16;
+    }
 
-        patch.x = (winSize.width  + block.x - 1) / block.x;
-        patch.y = (winSize.height + block.y - 1) / block.y;
+    patch.x = (winSize.width  + block.x - 1) / block.x;
+    patch.y = (winSize.height + block.y - 1) / block.y;
 
-        block.z = patch.z = 1;
-    }
+    block.z = patch.z = 1;
+}
 }
 
 inline int divUp(int total, int grain)
@@ -530,7 +531,7 @@ void arithmetic_run(const oclMat &src1, oclMat &dst, string kernelName, const ch
 
 void multiply_cus(const oclMat &src1, oclMat &dst, float scalar)
 {
-    arithmetic_run(src1, dst, "arithm_muls", &pyrlk, (void *)(&scalar));
+    arithmetic_run(src1, dst, "arithm_muls", &arithm_mul, (void *)(&scalar));
 }
 
 void pyrdown_run_cus(const oclMat &src, const oclMat &dst)
@@ -581,26 +582,26 @@ void pyrDown_cus(const oclMat &src, oclMat &dst)
 //
 //void callF(const oclMat& src, oclMat& dst, MultiplyScalar op, int mask)
 //{
-//  Mat srcTemp;
-//  Mat dstTemp;
-//  src.download(srcTemp);
-//  dst.download(dstTemp);
+//     Mat srcTemp;
+//     Mat dstTemp;
+//     src.download(srcTemp);
+//     dst.download(dstTemp);
 //
-//  int i;
-//  int j;
-//  int k;
-//  for(i = 0; i < srcTemp.rows; i++)
-//  {
-//      for(j = 0; j < srcTemp.cols; j++)
-//      {
-//          for(k = 0; k < srcTemp.channels(); k++)
-//          {
-//              ((float*)dstTemp.data)[srcTemp.channels() * (i * srcTemp.rows + j) + k] = (float)op(((float*)srcTemp.data)[srcTemp.channels() * (i * srcTemp.rows + j) + k]);
-//          }
-//      }
-//  }
+//     int i;
+//     int j;
+//     int k;
+//     for(i = 0; i < srcTemp.rows; i++)
+//     {
+//             for(j = 0; j < srcTemp.cols; j++)
+//             {
+//                     for(k = 0; k < srcTemp.channels(); k++)
+//                     {
+//                             ((float*)dstTemp.data)[srcTemp.channels() * (i * srcTemp.rows + j) + k] = (float)op(((float*)srcTemp.data)[srcTemp.channels() * (i * srcTemp.rows + j) + k]);
+//                     }
+//             }
+//     }
 //
-//  dst = dstTemp;
+//     dst = dstTemp;
 //}
 //
 //static inline bool isAligned(const unsigned char* ptr, size_t size)
@@ -622,54 +623,54 @@ void pyrDown_cus(const oclMat &src, oclMat &dst)
 //        return;
 //    }
 //
-//  Mat srcTemp;
-//  Mat dstTemp;
-//  src.download(srcTemp);
-//  dst.download(dstTemp);
+//     Mat srcTemp;
+//     Mat dstTemp;
+//     src.download(srcTemp);
+//     dst.download(dstTemp);
 //
-//  int x_shifted;
+//     int x_shifted;
 //
-//  int i;
-//  int j;
-//  for(i = 0; i < srcTemp.rows; i++)
-//  {
-//      const double* srcRow = (const double*)srcTemp.data + i * srcTemp.rows;
+//     int i;
+//     int j;
+//     for(i = 0; i < srcTemp.rows; i++)
+//     {
+//             const double* srcRow = (const double*)srcTemp.data + i * srcTemp.rows;
 //        double* dstRow = (double*)dstTemp.data + i * dstTemp.rows;;
 //
-//      for(j = 0; j < srcTemp.cols; j++)
-//      {
-//          x_shifted = j * 4;
+//             for(j = 0; j < srcTemp.cols; j++)
+//             {
+//                     x_shifted = j * 4;
 //
-//          if(x_shifted + 4 - 1 < srcTemp.cols)
-//          {
-//              dstRow[x_shifted    ] = op(srcRow[x_shifted    ]);
-//              dstRow[x_shifted + 1] = op(srcRow[x_shifted + 1]);
-//              dstRow[x_shifted + 2] = op(srcRow[x_shifted + 2]);
-//              dstRow[x_shifted + 3] = op(srcRow[x_shifted + 3]);
-//          }
-//          else
-//          {
-//              for (int real_x = x_shifted; real_x < srcTemp.cols; ++real_x)
-//              {
-//                  ((float*)dstTemp.data)[i * srcTemp.rows + real_x] = op(((float*)srcTemp.data)[i * srcTemp.rows + real_x]);
-//              }
-//          }
-//      }
-//  }
+//                     if(x_shifted + 4 - 1 < srcTemp.cols)
+//                     {
+//                             dstRow[x_shifted    ] = op(srcRow[x_shifted    ]);
+//                             dstRow[x_shifted + 1] = op(srcRow[x_shifted + 1]);
+//                             dstRow[x_shifted + 2] = op(srcRow[x_shifted + 2]);
+//                             dstRow[x_shifted + 3] = op(srcRow[x_shifted + 3]);
+//                     }
+//                     else
+//                     {
+//                             for (int real_x = x_shifted; real_x < srcTemp.cols; ++real_x)
+//                             {
+//                                     ((float*)dstTemp.data)[i * srcTemp.rows + real_x] = op(((float*)srcTemp.data)[i * srcTemp.rows + real_x]);
+//                             }
+//                     }
+//             }
+//     }
 //}
 //
 //void multiply(const oclMat& src1, double val, oclMat& dst, double scale = 1.0f);
 //void multiply(const oclMat& src1, double val, oclMat& dst, double scale)
 //{
 //    MultiplyScalar op(val, scale);
-//  //if(src1.channels() == 1 && dst.channels() == 1)
-//  //{
-//  //    callT(src1, dst, op, 0);
-//  //}
-//  //else
-//  //{
-//      callF(src1, dst, op, 0);
-//  //}
+//     //if(src1.channels() == 1 && dst.channels() == 1)
+//     //{
+//     //    callT(src1, dst, op, 0);
+//     //}
+//     //else
+//     //{
+//         callF(src1, dst, op, 0);
+//     //}
 //}
 
 cl_mem bindTexture(const oclMat &mat, int depth, int channels)
@@ -735,46 +736,69 @@ void releaseTexture(cl_mem texture)
 }
 
 void lkSparse_run(oclMat &I, oclMat &J,
-                  const oclMat &prevPts, oclMat &nextPts, oclMat &status, oclMat *err, bool GET_MIN_EIGENVALS, int ptcount,
+                  const oclMat &prevPts, oclMat &nextPts, oclMat &status, oclMat& err, bool /*GET_MIN_EIGENVALS*/, int ptcount,
                   int level, /*dim3 block, */dim3 patch, Size winSize, int iters)
 {
     Context  *clCxt = I.clCxt;
+    char platform[256] = {0};
+    cl_platform_id pid;
+    clGetDeviceInfo(*clCxt->impl->devices, CL_DEVICE_PLATFORM, sizeof(pid), &pid, NULL);
+    clGetPlatformInfo(pid, CL_PLATFORM_NAME, 256, platform, NULL);
+    std::string namestr = platform;
+    bool isImageSupported = true;
+    if(namestr.find("NVIDIA")!=string::npos || namestr.find("Intel")!=string::npos)
+        isImageSupported = false;
+
+    int elemCntPerRow = I.step / I.elemSize();
 
     string kernelName = "lkSparse";
 
-    size_t localThreads[3]  = { 8, 32, 1 };
-    size_t globalThreads[3] = { 8 * ptcount, 32, 1};
+
+    size_t localThreads[3]  = { 8, isImageSupported?8:32, 1 };
+    size_t globalThreads[3] = { 8 * ptcount, isImageSupported?8:32, 1};
 
     int cn = I.oclchannels();
 
-    bool calcErr;
-    if (err)
+    char calcErr;
+    if (level == 0)
     {
-        calcErr = true;
+        calcErr = 1;
     }
     else
     {
-        calcErr = false;
+        calcErr = 0;
     }
-    calcErr = true;
-
-    cl_mem ITex = bindTexture(I, I.depth(), cn);
-    cl_mem JTex = bindTexture(J, J.depth(), cn);
 
     vector<pair<size_t , const void *> > args;
+    cl_mem ITex;
+    cl_mem JTex;
+    if (isImageSupported)
+    {
+        ITex = bindTexture(I, I.depth(), cn);
+        JTex = bindTexture(J, J.depth(), cn);
+    }
+    else
+    {
+        ITex = (cl_mem)I.data;
+        JTex = (cl_mem)J.data;
+    }
 
     args.push_back( make_pair( sizeof(cl_mem), (void *)&ITex ));
     args.push_back( make_pair( sizeof(cl_mem), (void *)&JTex ));
-
+    //cl_mem clmD = clCreateBuffer(clCxt, CL_MEM_READ_WRITE, ptcount * sizeof(float), NULL, NULL);
     args.push_back( make_pair( sizeof(cl_mem), (void *)&prevPts.data ));
     args.push_back( make_pair( sizeof(cl_int), (void *)&prevPts.step ));
     args.push_back( make_pair( sizeof(cl_mem), (void *)&nextPts.data ));
     args.push_back( make_pair( sizeof(cl_int), (void *)&nextPts.step ));
     args.push_back( make_pair( sizeof(cl_mem), (void *)&status.data ));
-    //args.push_back( make_pair( sizeof(cl_mem), (void *)&(err->data) ));
+    args.push_back( make_pair( sizeof(cl_mem), (void *)&err.data ));
     args.push_back( make_pair( sizeof(cl_int), (void *)&level ));
     args.push_back( make_pair( sizeof(cl_int), (void *)&I.rows ));
     args.push_back( make_pair( sizeof(cl_int), (void *)&I.cols ));
+    if (!isImageSupported)
+    {
+        args.push_back( make_pair( sizeof(cl_int), (void *)&elemCntPerRow ) );
+    }
     args.push_back( make_pair( sizeof(cl_int), (void *)&patch.x ));
     args.push_back( make_pair( sizeof(cl_int), (void *)&patch.y ));
     args.push_back( make_pair( sizeof(cl_int), (void *)&cn ));
@@ -782,27 +806,29 @@ void lkSparse_run(oclMat &I, oclMat &J,
     args.push_back( make_pair( sizeof(cl_int), (void *)&winSize.height ));
     args.push_back( make_pair( sizeof(cl_int), (void *)&iters ));
     args.push_back( make_pair( sizeof(cl_char), (void *)&calcErr ));
-    args.push_back( make_pair( sizeof(cl_char), (void *)&GET_MIN_EIGENVALS ));
+    //args.push_back( make_pair( sizeof(cl_char), (void *)&GET_MIN_EIGENVALS ));
 
-    openCLExecuteKernel2(clCxt, &pyrlk, kernelName, globalThreads, localThreads, args, I.oclchannels(), I.depth(), CLFLUSH);
+    if (isImageSupported)
+    {
+        openCLExecuteKernel2(clCxt, &pyrlk, kernelName, globalThreads, localThreads, args, I.oclchannels(), I.depth(), CLFLUSH);
 
-    releaseTexture(ITex);
-    releaseTexture(JTex);
+        releaseTexture(ITex);
+        releaseTexture(JTex);
+    }
+    else
+    {
+        //printf("Warning: The image2d_t is not supported by the device. Using alternative method!\n");
+        openCLExecuteKernel2(clCxt, &pyrlk_no_image, kernelName, globalThreads, localThreads, args, I.oclchannels(), I.depth(), CLFLUSH);
+    }
 }
 
 void cv::ocl::PyrLKOpticalFlow::sparse(const oclMat &prevImg, const oclMat &nextImg, const oclMat &prevPts, oclMat &nextPts, oclMat &status, oclMat *err)
 {
-    if (prevImg.clCxt->impl->devName.find("Intel(R) HD Graphics") != string::npos)
-    {
-        cout << " Intel HD GPU device unsupported " << endl;
-        return;
-    }
-
     if (prevPts.empty())
     {
         nextPts.release();
         status.release();
-        if (err) err->release();
+        //if (err) err->release();
         return;
     }
 
@@ -836,8 +862,15 @@ void cv::ocl::PyrLKOpticalFlow::sparse(const oclMat &prevImg, const oclMat &next
     //status.setTo(Scalar::all(1));
     setTo(status, Scalar::all(1));
 
-    //if (err)
-    //    ensureSizeIsEnough(1, prevPts.cols, CV_32FC1, *err);
+    bool errMat = false;
+    if (!err)
+    {
+        err = new oclMat(1, prevPts.cols, CV_32FC1);
+        errMat = true;
+    }
+    else
+        ensureSizeIsEnough(1, prevPts.cols, CV_32FC1, *err);
+    //ensureSizeIsEnough(1, prevPts.cols, CV_32FC1, err);
 
     // build the image pyramids.
 
@@ -872,17 +905,22 @@ void cv::ocl::PyrLKOpticalFlow::sparse(const oclMat &prevImg, const oclMat &next
     for (int level = maxLevel; level >= 0; level--)
     {
         lkSparse_run(prevPyr_[level], nextPyr_[level],
-                     prevPts, nextPts, status, level == 0 && err ? err : 0, getMinEigenVals, prevPts.cols,
+                     prevPts, nextPts, status, *err, getMinEigenVals, prevPts.cols,
                      level, /*block, */patch, winSize, iters);
     }
 
     clFinish(prevImg.clCxt->impl->clCmdQueue);
+
+    if(errMat)
+        delete err;
 }
 
 void lkDense_run(oclMat &I, oclMat &J, oclMat &u, oclMat &v,
                  oclMat &prevU, oclMat &prevV, oclMat *err, Size winSize, int iters)
 {
     Context  *clCxt = I.clCxt;
+    bool isImageSupported = clCxt->impl->devName.find("Intel(R) HD Graphics") == string::npos;
+    int elemCntPerRow = I.step / I.elemSize();
 
     string kernelName = "lkDense";
 
@@ -901,8 +939,19 @@ void lkDense_run(oclMat &I, oclMat &J, oclMat &u, oclMat &v,
         calcErr = false;
     }
 
-    cl_mem ITex = bindTexture(I, I.depth(), cn);
-    cl_mem JTex = bindTexture(J, J.depth(), cn);
+    cl_mem ITex;
+    cl_mem JTex;
+
+    if (isImageSupported)
+    {
+        ITex = bindTexture(I, I.depth(), cn);
+        JTex = bindTexture(J, J.depth(), cn);
+    }
+    else
+    {
+        ITex = (cl_mem)I.data;
+        JTex = (cl_mem)J.data;
+    }
 
     //int2 halfWin = {(winSize.width - 1) / 2, (winSize.height - 1) / 2};
     //const int patchWidth  = 16 + 2 * halfWin.x;
@@ -926,15 +975,27 @@ void lkDense_run(oclMat &I, oclMat &J, oclMat &u, oclMat &v,
     args.push_back( make_pair( sizeof(cl_int), (void *)&I.cols ));
     //args.push_back( make_pair( sizeof(cl_mem), (void *)&(*err).data ));
     //args.push_back( make_pair( sizeof(cl_int), (void *)&(*err).step ));
+    if (!isImageSupported)
+    {
+        args.push_back( make_pair( sizeof(cl_int), (void *)&elemCntPerRow ) );
+    }
     args.push_back( make_pair( sizeof(cl_int), (void *)&winSize.width ));
     args.push_back( make_pair( sizeof(cl_int), (void *)&winSize.height ));
     args.push_back( make_pair( sizeof(cl_int), (void *)&iters ));
     args.push_back( make_pair( sizeof(cl_char), (void *)&calcErr ));
 
-    openCLExecuteKernel2(clCxt, &pyrlk, kernelName, globalThreads, localThreads, args, I.oclchannels(), I.depth(), CLFLUSH);
+    if (isImageSupported)
+    {
+        openCLExecuteKernel2(clCxt, &pyrlk, kernelName, globalThreads, localThreads, args, I.oclchannels(), I.depth(), CLFLUSH);
 
-    releaseTexture(ITex);
-    releaseTexture(JTex);
+        releaseTexture(ITex);
+        releaseTexture(JTex);
+    }
+    else
+    {
+        //printf("Warning: The image2d_t is not supported by the device. Using alternative method!\n");
+        openCLExecuteKernel2(clCxt, &pyrlk_no_image, kernelName, globalThreads, localThreads, args, I.oclchannels(), I.depth(), CLFLUSH);
+    }
 }
 
 void cv::ocl::PyrLKOpticalFlow::dense(const oclMat &prevImg, const oclMat &nextImg, oclMat &u, oclMat &v, oclMat *err)
index b594a34..7c747ee 100644 (file)
@@ -118,9 +118,9 @@ TEST_P(Sparse, Mat)
     cv::Mat status_mat(1, d_status.cols, CV_8UC1, (void *)&status[0]);
     d_status.download(status_mat);
 
-    //std::vector<float> err(d_err.cols);
-    //cv::Mat err_mat(1, d_err.cols, CV_32FC1, (void*)&err[0]);
-    //d_err.download(err_mat);
+    std::vector<float> err(d_err.cols);
+    cv::Mat err_mat(1, d_err.cols, CV_32FC1, (void*)&err[0]);
+    d_err.download(err_mat);
 
     std::vector<cv::Point2f> nextPts_gold;
     std::vector<unsigned char> status_gold;
@@ -153,9 +153,9 @@ TEST_P(Sparse, Mat)
         }
     }
 
-    double bad_ratio = static_cast<double>(mistmatch) / (nextPts.size() * 2);
+    double bad_ratio = static_cast<double>(mistmatch) / (nextPts.size());
 
-    ASSERT_LE(bad_ratio, 0.05f);
+    ASSERT_LE(bad_ratio, 0.02f);
 
 }