Merge pull request #2063 from vpisarev:ocl_moments

modules/imgproc/src/moments.cpp

@@ -39,6 +39,7 @@
 //
 //M*/
 #include "precomp.hpp"
+#include "opencl_kernels.hpp"
 
 namespace cv
 {
@@ -362,106 +363,175 @@ Moments::Moments( double _m00, double _m10, double _m01, double _m20, double _m1
     nu30 = mu30*s3; nu21 = mu21*s3; nu12 = mu12*s3; nu03 = mu03*s3;
 }
 
+static bool ocl_moments( InputArray _src, Moments& m)
+{
+    const int TILE_SIZE = 32;
+    const int K = 10;
+    ocl::Kernel k("moments", ocl::imgproc::moments_oclsrc, format("-D TILE_SIZE=%d", TILE_SIZE));
+    if( k.empty() )
+        return false;
+
+    UMat src = _src.getUMat();
+    Size sz = src.size();
+    int xtiles = (sz.width + TILE_SIZE-1)/TILE_SIZE;
+    int ytiles = (sz.height + TILE_SIZE-1)/TILE_SIZE;
+    int ntiles = xtiles*ytiles;
+    UMat umbuf(1, ntiles*K, CV_32S);
+
+    size_t globalsize[] = {xtiles, sz.height}, localsize[] = {1, TILE_SIZE};
+    bool ok = k.args(ocl::KernelArg::ReadOnly(src),
+                     ocl::KernelArg::PtrWriteOnly(umbuf),
+                     xtiles).run(2, globalsize, localsize, true);
+    if(!ok)
+        return false;
+    Mat mbuf = umbuf.getMat(ACCESS_READ);
+    for( int i = 0; i < ntiles; i++ )
+    {
+        double x = (i % xtiles)*TILE_SIZE, y = (i / xtiles)*TILE_SIZE;
+        const int* mom = mbuf.ptr<int>() + i*K;
+        double xm = x * mom[0], ym = y * mom[0];
+
+        // accumulate moments computed in each tile
+
+        // + m00 ( = m00' )
+        m.m00 += mom[0];
+
+        // + m10 ( = m10' + x*m00' )
+        m.m10 += mom[1] + xm;
+
+        // + m01 ( = m01' + y*m00' )
+        m.m01 += mom[2] + ym;
+
+        // + m20 ( = m20' + 2*x*m10' + x*x*m00' )
+        m.m20 += mom[3] + x * (mom[1] * 2 + xm);
+
+        // + m11 ( = m11' + x*m01' + y*m10' + x*y*m00' )
+        m.m11 += mom[4] + x * (mom[2] + ym) + y * mom[1];
+
+        // + m02 ( = m02' + 2*y*m01' + y*y*m00' )
+        m.m02 += mom[5] + y * (mom[2] * 2 + ym);
+
+        // + m30 ( = m30' + 3*x*m20' + 3*x*x*m10' + x*x*x*m00' )
+        m.m30 += mom[6] + x * (3. * mom[3] + x * (3. * mom[1] + xm));
+
+        // + m21 ( = m21' + x*(2*m11' + 2*y*m10' + x*m01' + x*y*m00') + y*m20')
+        m.m21 += mom[7] + x * (2 * (mom[4] + y * mom[1]) + x * (mom[2] + ym)) + y * mom[3];
+
+        // + m12 ( = m12' + y*(2*m11' + 2*x*m01' + y*m10' + x*y*m00') + x*m02')
+        m.m12 += mom[8] + y * (2 * (mom[4] + x * mom[2]) + y * (mom[1] + xm)) + x * mom[5];
+
+        // + m03 ( = m03' + 3*y*m02' + 3*y*y*m01' + y*y*y*m00' )
+        m.m03 += mom[9] + y * (3. * mom[5] + y * (3. * mom[2] + ym));
+    }
+
+    return true;
+}
+
 }
 
 
 cv::Moments cv::moments( InputArray _src, bool binary )
 {
     const int TILE_SIZE = 32;
-    Mat mat = _src.getMat();
     MomentsInTileFunc func = 0;
     uchar nzbuf[TILE_SIZE*TILE_SIZE];
     Moments m;
-    int type = mat.type();
+    int type = _src.type();
     int depth = CV_MAT_DEPTH( type );
     int cn = CV_MAT_CN( type );
-
-    if( mat.checkVector(2) >= 0 && (depth == CV_32F || depth == CV_32S))
-        return contourMoments(mat);
-
-    Size size = mat.size();
-
-    if( cn > 1 )
-        CV_Error( CV_StsBadArg, "Invalid image type" );
+    Size size = _src.size();
 
     if( size.width <= 0 || size.height <= 0 )
         return m;
 
-    if( binary || depth == CV_8U )
-        func = momentsInTile<uchar, int, int>;
-    else if( depth == CV_16U )
-        func = momentsInTile<ushort, int, int64>;
-    else if( depth == CV_16S )
-        func = momentsInTile<short, int, int64>;
-    else if( depth == CV_32F )
-        func = momentsInTile<float, double, double>;
-    else if( depth == CV_64F )
-        func = momentsInTile<double, double, double>;
+    if( ocl::useOpenCL() && type == CV_8UC1 && !binary &&
+        _src.isUMat() && ocl_moments(_src, m) )
+        ;
     else
-        CV_Error( CV_StsUnsupportedFormat, "" );
-
-    Mat src0(mat);
-
-    for( int y = 0; y < size.height; y += TILE_SIZE )
     {
-        Size tileSize;
-        tileSize.height = std::min(TILE_SIZE, size.height - y);
+        Mat mat = _src.getMat();
+        if( mat.checkVector(2) >= 0 && (depth == CV_32F || depth == CV_32S))
+            return contourMoments(mat);
+
+        if( cn > 1 )
+            CV_Error( CV_StsBadArg, "Invalid image type (must be single-channel)" );
+
+        if( binary || depth == CV_8U )
+            func = momentsInTile<uchar, int, int>;
+        else if( depth == CV_16U )
+            func = momentsInTile<ushort, int, int64>;
+        else if( depth == CV_16S )
+            func = momentsInTile<short, int, int64>;
+        else if( depth == CV_32F )
+            func = momentsInTile<float, double, double>;
+        else if( depth == CV_64F )
+            func = momentsInTile<double, double, double>;
+        else
+            CV_Error( CV_StsUnsupportedFormat, "" );
+
+        Mat src0(mat);
 
-        for( int x = 0; x < size.width; x += TILE_SIZE )
+        for( int y = 0; y < size.height; y += TILE_SIZE )
         {
-            tileSize.width = std::min(TILE_SIZE, size.width - x);
-            Mat src(src0, cv::Rect(x, y, tileSize.width, tileSize.height));
+            Size tileSize;
+            tileSize.height = std::min(TILE_SIZE, size.height - y);
 
-            if( binary )
+            for( int x = 0; x < size.width; x += TILE_SIZE )
             {
-                cv::Mat tmp(tileSize, CV_8U, nzbuf);
-                cv::compare( src, 0, tmp, CV_CMP_NE );
-                src = tmp;
-            }
+                tileSize.width = std::min(TILE_SIZE, size.width - x);
+                Mat src(src0, cv::Rect(x, y, tileSize.width, tileSize.height));
 
-            double mom[10];
-            func( src, mom );
+                if( binary )
+                {
+                    cv::Mat tmp(tileSize, CV_8U, nzbuf);
+                    cv::compare( src, 0, tmp, CV_CMP_NE );
+                    src = tmp;
+                }
 
-            if(binary)
-            {
-                double s = 1./255;
-                for( int k = 0; k < 10; k++ )
-                    mom[k] *= s;
-            }
+                double mom[10];
+                func( src, mom );
+
+                if(binary)
+                {
+                    double s = 1./255;
+                    for( int k = 0; k < 10; k++ )
+                        mom[k] *= s;
+                }
 
-            double xm = x * mom[0], ym = y * mom[0];
+                double xm = x * mom[0], ym = y * mom[0];
 
-            // accumulate moments computed in each tile
+                // accumulate moments computed in each tile
 
-            // + m00 ( = m00' )
-            m.m00 += mom[0];
+                // + m00 ( = m00' )
+                m.m00 += mom[0];
 
-            // + m10 ( = m10' + x*m00' )
-            m.m10 += mom[1] + xm;
+                // + m10 ( = m10' + x*m00' )
+                m.m10 += mom[1] + xm;
 
-            // + m01 ( = m01' + y*m00' )
-            m.m01 += mom[2] + ym;
+                // + m01 ( = m01' + y*m00' )
+                m.m01 += mom[2] + ym;
 
-            // + m20 ( = m20' + 2*x*m10' + x*x*m00' )
-            m.m20 += mom[3] + x * (mom[1] * 2 + xm);
+                // + m20 ( = m20' + 2*x*m10' + x*x*m00' )
+                m.m20 += mom[3] + x * (mom[1] * 2 + xm);
 
-            // + m11 ( = m11' + x*m01' + y*m10' + x*y*m00' )
-            m.m11 += mom[4] + x * (mom[2] + ym) + y * mom[1];
+                // + m11 ( = m11' + x*m01' + y*m10' + x*y*m00' )
+                m.m11 += mom[4] + x * (mom[2] + ym) + y * mom[1];
 
-            // + m02 ( = m02' + 2*y*m01' + y*y*m00' )
-            m.m02 += mom[5] + y * (mom[2] * 2 + ym);
+                // + m02 ( = m02' + 2*y*m01' + y*y*m00' )
+                m.m02 += mom[5] + y * (mom[2] * 2 + ym);
 
-            // + m30 ( = m30' + 3*x*m20' + 3*x*x*m10' + x*x*x*m00' )
-            m.m30 += mom[6] + x * (3. * mom[3] + x * (3. * mom[1] + xm));
+                // + m30 ( = m30' + 3*x*m20' + 3*x*x*m10' + x*x*x*m00' )
+                m.m30 += mom[6] + x * (3. * mom[3] + x * (3. * mom[1] + xm));
 
-            // + m21 ( = m21' + x*(2*m11' + 2*y*m10' + x*m01' + x*y*m00') + y*m20')
-            m.m21 += mom[7] + x * (2 * (mom[4] + y * mom[1]) + x * (mom[2] + ym)) + y * mom[3];
+                // + m21 ( = m21' + x*(2*m11' + 2*y*m10' + x*m01' + x*y*m00') + y*m20')
+                m.m21 += mom[7] + x * (2 * (mom[4] + y * mom[1]) + x * (mom[2] + ym)) + y * mom[3];
 
-            // + m12 ( = m12' + y*(2*m11' + 2*x*m01' + y*m10' + x*y*m00') + x*m02')
-            m.m12 += mom[8] + y * (2 * (mom[4] + x * mom[2]) + y * (mom[1] + xm)) + x * mom[5];
+                // + m12 ( = m12' + y*(2*m11' + 2*x*m01' + y*m10' + x*y*m00') + x*m02')
+                m.m12 += mom[8] + y * (2 * (mom[4] + x * mom[2]) + y * (mom[1] + xm)) + x * mom[5];
 
-            // + m03 ( = m03' + 3*y*m02' + 3*y*y*m01' + y*y*y*m00' )
-            m.m03 += mom[9] + y * (3. * mom[5] + y * (3. * mom[2] + ym));
+                // + m03 ( = m03' + 3*y*m02' + 3*y*y*m01' + y*y*y*m00' )
+                m.m03 += mom[9] + y * (3. * mom[5] + y * (3. * mom[2] + ym));
+            }
         }
     }
 

modules/imgproc/src/opencl/moments.cl

+/* See LICENSE file in the root OpenCV directory */
+
+#if TILE_SIZE != 32
+#error "TILE SIZE should be 32"
+#endif
+
+__kernel void moments(__global const uchar* src, int src_step, int src_offset,
+                      int src_rows, int src_cols, __global int* mom0, int xtiles)
+{
+    int x0 = get_global_id(0);
+    int y0 = get_group_id(1);
+    int x, y = get_local_id(1);
+    int x_min = x0*TILE_SIZE;
+    int ypix = y0*TILE_SIZE + y;
+    __local int mom[TILE_SIZE][10];
+
+    if( x_min < src_cols && y0*TILE_SIZE < src_rows )
+    {
+        if( ypix < src_rows )
+        {
+            int x_max = min(src_cols - x_min, TILE_SIZE);
+            __global const uchar* ptr = src + src_offset + ypix*src_step + x_min;
+            int4 S = (int4)(0,0,0,0), p;
+
+            #define SUM_ELEM(elem, ofs) \
+                (int4)(1, (ofs), (ofs)*(ofs), (ofs)*(ofs)*(ofs))*elem
+
+            x = x_max & -4;
+            if( x_max >= 4 )
+            {
+                p = convert_int4(vload4(0, ptr));
+                S += SUM_ELEM(p.s0, 0) + SUM_ELEM(p.s1, 1) + SUM_ELEM(p.s2, 2) + SUM_ELEM(p.s3, 3);
+
+                if( x_max >= 8 )
+                {
+                    p = convert_int4(vload4(0, ptr+4));
+                    S += SUM_ELEM(p.s0, 4) + SUM_ELEM(p.s1, 5) + SUM_ELEM(p.s2, 6) + SUM_ELEM(p.s3, 7);
+
+                    if( x_max >= 12 )
+                    {
+                        p = convert_int4(vload4(0, ptr+8));
+                        S += SUM_ELEM(p.s0, 8) + SUM_ELEM(p.s1, 9) + SUM_ELEM(p.s2, 10) + SUM_ELEM(p.s3, 11);
+
+                        if( x_max >= 16 )
+                        {
+                            p = convert_int4(vload4(0, ptr+12));
+                            S += SUM_ELEM(p.s0, 12) + SUM_ELEM(p.s1, 13) + SUM_ELEM(p.s2, 14) + SUM_ELEM(p.s3, 15);
+                        }
+                    }
+                }
+            }
+
+            if( x_max >= 20 )
+            {
+                p = convert_int4(vload4(0, ptr+16));
+                S += SUM_ELEM(p.s0, 16) + SUM_ELEM(p.s1, 17) + SUM_ELEM(p.s2, 18) + SUM_ELEM(p.s3, 19);
+
+                if( x_max >= 24 )
+                {
+                    p = convert_int4(vload4(0, ptr+20));
+                    S += SUM_ELEM(p.s0, 20) + SUM_ELEM(p.s1, 21) + SUM_ELEM(p.s2, 22) + SUM_ELEM(p.s3, 23);
+
+                    if( x_max >= 28 )
+                    {
+                        p = convert_int4(vload4(0, ptr+24));
+                        S += SUM_ELEM(p.s0, 24) + SUM_ELEM(p.s1, 25) + SUM_ELEM(p.s2, 26) + SUM_ELEM(p.s3, 27);
+
+                        if( x_max >= 32 )
+                        {
+                            p = convert_int4(vload4(0, ptr+28));
+                            S += SUM_ELEM(p.s0, 28) + SUM_ELEM(p.s1, 29) + SUM_ELEM(p.s2, 30) + SUM_ELEM(p.s3, 31);
+                        }
+                    }
+                }
+            }
+
+            if( x < x_max )
+            {
+                int ps = ptr[x];
+                S += SUM_ELEM(ps, x);
+                if( x+1 < x_max )
+                {
+                    ps = ptr[x+1];
+                    S += SUM_ELEM(ps, x+1);
+                    if( x+2 < x_max )
+                    {
+                        ps = ptr[x+2];
+                        S += SUM_ELEM(ps, x+2);
+                    }
+                }
+            }
+
+            int sy = y*y;
+
+            mom[y][0] = S.s0;
+            mom[y][1] = S.s1;
+            mom[y][2] = y*S.s0;
+            mom[y][3] = S.s2;
+            mom[y][4] = y*S.s1;
+            mom[y][5] = sy*S.s0;
+            mom[y][6] = S.s3;
+            mom[y][7] = y*S.s2;
+            mom[y][8] = sy*S.s1;
+            mom[y][9] = y*sy*S.s0;
+        }
+        else
+            mom[y][0] = mom[y][1] = mom[y][2] = mom[y][3] = mom[y][4] =
+                mom[y][5] = mom[y][6] = mom[y][7] = mom[y][8] = mom[y][9] = 0;
+        barrier(CLK_LOCAL_MEM_FENCE);
+
+        #define REDUCE(d) \
+        if( y < d ) \
+        { \
+            mom[y][0] += mom[y+d][0]; \
+            mom[y][1] += mom[y+d][1]; \
+            mom[y][2] += mom[y+d][2]; \
+            mom[y][3] += mom[y+d][3]; \
+            mom[y][4] += mom[y+d][4]; \
+            mom[y][5] += mom[y+d][5]; \
+            mom[y][6] += mom[y+d][6]; \
+            mom[y][7] += mom[y+d][7]; \
+            mom[y][8] += mom[y+d][8]; \
+            mom[y][9] += mom[y+d][9]; \
+        } \
+        barrier(CLK_LOCAL_MEM_FENCE)
+
+        REDUCE(16);
+        REDUCE(8);
+        REDUCE(4);
+        REDUCE(2);
+
+        if( y == 0 )
+        {
+            __global int* momout = mom0 + (y0*xtiles + x0)*10;
+            momout[0] = mom[0][0] + mom[1][0];
+            momout[1] = mom[0][1] + mom[1][1];
+            momout[2] = mom[0][2] + mom[1][2];
+            momout[3] = mom[0][3] + mom[1][3];
+            momout[4] = mom[0][4] + mom[1][4];
+            momout[5] = mom[0][5] + mom[1][5];
+            momout[6] = mom[0][6] + mom[1][6];
+            momout[7] = mom[0][7] + mom[1][7];
+            momout[8] = mom[0][8] + mom[1][8];
+            momout[9] = mom[0][9] + mom[1][9];
+        }
+    }
+}

modules/imgproc/test/test_moments.cpp

@@ -43,6 +43,13 @@
 using namespace cv;
 using namespace std;
 
+#define OCL_TUNING_MODE 0
+#if OCL_TUNING_MODE
+#define OCL_TUNING_MODE_ONLY(code) code
+#else
+#define OCL_TUNING_MODE_ONLY(code)
+#endif
+
 // image moments
 class CV_MomentsTest : public cvtest::ArrayTest
 {
@@ -60,6 +67,7 @@ protected:
     void run_func();
     int coi;
     bool is_binary;
+    bool try_umat;
 };
 
 
@@ -70,6 +78,7 @@ CV_MomentsTest::CV_MomentsTest()
     test_array[REF_OUTPUT].push_back(NULL);
     coi = -1;
     is_binary = false;
+    OCL_TUNING_MODE_ONLY(test_case_count = 10);
     //element_wise_relative_error = false;
 }
 
@@ -96,25 +105,38 @@ void CV_MomentsTest::get_minmax_bounds( int i, int j, int type, Scalar& low, Sca
     }
 }
 
-
 void CV_MomentsTest::get_test_array_types_and_sizes( int test_case_idx,
                                                 vector<vector<Size> >& sizes, vector<vector<int> >& types )
 {
     RNG& rng = ts->get_rng();
     cvtest::ArrayTest::get_test_array_types_and_sizes( test_case_idx, sizes, types );
-    int cn = cvtest::randInt(rng) % 4 + 1;
+    int cn = (cvtest::randInt(rng) % 4) + 1;
     int depth = cvtest::randInt(rng) % 4;
     depth = depth == 0 ? CV_8U : depth == 1 ? CV_16U : depth == 2 ? CV_16S : CV_32F;
-    if( cn == 2 )
+
+    is_binary = cvtest::randInt(rng) % 2 != 0;
+    if( depth == 0 && !is_binary )
+        try_umat = cvtest::randInt(rng) % 5 != 0;
+    else
+        try_umat = cvtest::randInt(rng) % 2 != 0;
+
+    if( cn == 2 || try_umat )
         cn = 1;
 
+    OCL_TUNING_MODE_ONLY(
+    cn = 1;
+    depth = CV_8U;
+    try_umat = true;
+    is_binary = false;
+    sizes[INPUT][0] = Size(1024,768)
+    );
+
     types[INPUT][0] = CV_MAKETYPE(depth, cn);
     types[OUTPUT][0] = types[REF_OUTPUT][0] = CV_64FC1;
     sizes[OUTPUT][0] = sizes[REF_OUTPUT][0] = cvSize(MOMENT_COUNT,1);
     if(CV_MAT_DEPTH(types[INPUT][0])>=CV_32S)
         sizes[INPUT][0].width = MAX(sizes[INPUT][0].width, 3);
 
-    is_binary = cvtest::randInt(rng) % 2 != 0;
     coi = 0;
     cvmat_allowed = true;
     if( cn > 1 )
@@ -149,7 +171,25 @@ void CV_MomentsTest::run_func()
 {
     CvMoments* m = (CvMoments*)test_mat[OUTPUT][0].ptr<double>();
     double* others = (double*)(m + 1);
-    cvMoments( test_array[INPUT][0], m, is_binary );
+    if( try_umat )
+    {
+        UMat u;
+        test_mat[INPUT][0].clone().copyTo(u);
+        OCL_TUNING_MODE_ONLY(
+            static double ttime = 0;
+            static int ncalls = 0;
+            moments(u, is_binary != 0);
+            double t = (double)getTickCount());
+        Moments new_m = moments(u, is_binary != 0);
+        OCL_TUNING_MODE_ONLY(
+            ttime += (double)getTickCount() - t;
+            ncalls++;
+            printf("%g\n", ttime/ncalls/u.total()));
+        *m = new_m;
+    }
+    else
+        cvMoments( test_array[INPUT][0], m, is_binary );
+
     others[0] = cvGetNormalizedCentralMoment( m, 2, 0 );
     others[1] = cvGetNormalizedCentralMoment( m, 1, 1 );
     others[2] = cvGetNormalizedCentralMoment( m, 0, 2 );