modules/ocl/src/moments.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                           License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2010-2012, Institute Of Software Chinese Academy Of Science, all rights reserved.
  14 // Copyright (C) 2010-2012, Advanced Micro Devices, Inc., all rights reserved.
  15 // Copyright (C) 2010-2012, Multicoreware, Inc., all rights reserved.
  16 // Third party copyrights are property of their respective owners.
  17 //
  18 // @Authors
  19 //    Sen Liu, swjtuls1987@126.com
  20 //
  21 // Redistribution and use in source and binary forms, with or without modification,
  22 // are permitted provided that the following conditions are met:
  23 //
  24 //   * Redistribution's of source code must retain the above copyright notice,
  25 //     this list of conditions and the following disclaimer.
  26 //
  27 //   * Redistribution's in binary form must reproduce the above copyright notice,
  28 //     this list of conditions and the following disclaimer in the documentation
  29 //     and/or other oclMaterials provided with the distribution.
  30 //
  31 //   * The name of the copyright holders may not be used to endorse or promote products
  32 //     derived from this software without specific prior written permission.
  33 //
  34 // This software is provided by the copyright holders and contributors "as is" and
  35 // any express or implied warranties, including, but not limited to, the implied
  36 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  37 // In no event shall the Intel Corporation or contributors be liable for any direct,
  38 // indirect, incidental, special, exemplary, or consequential damages
  39 // (including, but not limited to, procurement of substitute goods or services;
  40 // loss of use, data, or profits; or business interruption) however caused
  41 // and on any theory of liability, whether in contract, strict liability,
  42 // or tort (including negligence or otherwise) arising in any way out of
  43 // the use of this software, even if advised of the possibility of such damage.
  44 //
  45 //M*/
  46 #include "precomp.hpp"
  47 #include <iostream>
  48 namespace cv
  49 {
  50 namespace ocl
  51 {
  52 extern const char *moments;
  53
  54 // The function calculates center of gravity and the central second order moments
  55 static void icvCompleteMomentState( CvMoments* moments )
  56 {
  57     double cx = 0, cy = 0;
  58     double mu20, mu11, mu02;
  59
  60     assert( moments != 0 );
  61     moments->inv_sqrt_m00 = 0;
  62
  63     if( fabs(moments->m00) > DBL_EPSILON )
  64     {
  65         double inv_m00 = 1. / moments->m00;
  66         cx = moments->m10 * inv_m00;
  67         cy = moments->m01 * inv_m00;
  68         moments->inv_sqrt_m00 = std::sqrt( fabs(inv_m00) );
  69     }
  70
  71     // mu20 = m20 - m10*cx
  72     mu20 = moments->m20 - moments->m10 * cx;
  73     // mu11 = m11 - m10*cy
  74     mu11 = moments->m11 - moments->m10 * cy;
  75     // mu02 = m02 - m01*cy
  76     mu02 = moments->m02 - moments->m01 * cy;
  77
  78     moments->mu20 = mu20;
  79     moments->mu11 = mu11;
  80     moments->mu02 = mu02;
  81
  82     // mu30 = m30 - cx*(3*mu20 + cx*m10)
  83     moments->mu30 = moments->m30 - cx * (3 * mu20 + cx * moments->m10);
  84     mu11 += mu11;
  85     // mu21 = m21 - cx*(2*mu11 + cx*m01) - cy*mu20
  86     moments->mu21 = moments->m21 - cx * (mu11 + cx * moments->m01) - cy * mu20;
  87     // mu12 = m12 - cy*(2*mu11 + cy*m10) - cx*mu02
  88     moments->mu12 = moments->m12 - cy * (mu11 + cy * moments->m10) - cx * mu02;
  89     // mu03 = m03 - cy*(3*mu02 + cy*m01)
  90     moments->mu03 = moments->m03 - cy * (3 * mu02 + cy * moments->m01);
  91 }
  92
  93
  94 static void icvContourMoments( CvSeq* contour, CvMoments* mom )
  95 {
  96     if( contour->total )
  97     {
  98         CvSeqReader reader;
  99         int lpt = contour->total;
 100         double a00, a10, a01, a20, a11, a02, a30, a21, a12, a03;
 101
 102         cvStartReadSeq( contour, &reader, 0 );
 103
 104         size_t reader_size = lpt << 1;
 105         cv::Mat reader_mat(1,reader_size,CV_32FC1);
 106
 107         bool is_float = CV_SEQ_ELTYPE(contour) == CV_32FC2;
 108
 109         if (!cv::ocl::Context::getContext()->supportsFeature(Context::CL_DOUBLE) && is_float)
 110         {
 111             CV_Error(CV_StsUnsupportedFormat, "Moments - double is not supported by your GPU!");
 112         }
 113
 114         if( is_float )
 115         {
 116             for(size_t i = 0; i < reader_size; ++i)
 117             {
 118                 reader_mat.at<float>(0, i++) = ((CvPoint2D32f*)(reader.ptr))->x;
 119                 reader_mat.at<float>(0, i) = ((CvPoint2D32f*)(reader.ptr))->y;
 120                 CV_NEXT_SEQ_ELEM( contour->elem_size, reader );
 121             }
 122         }
 123         else
 124         {
 125             for(size_t i = 0; i < reader_size; ++i)
 126             {
 127                 reader_mat.at<float>(0, i++) = ((CvPoint*)(reader.ptr))->x;
 128                 reader_mat.at<float>(0, i) = ((CvPoint*)(reader.ptr))->y;
 129                 CV_NEXT_SEQ_ELEM( contour->elem_size, reader );
 130             }
 131         }
 132
 133         cv::ocl::oclMat dst_a(10, lpt, CV_64FC1);
 134         cv::ocl::oclMat reader_oclmat(reader_mat);
 135         int llength = std::min(lpt,128);
 136         size_t localThreads[3]  = { llength, 1, 1};
 137         size_t globalThreads[3] = { lpt, 1, 1};
 138         vector<pair<size_t , const void *> > args;
 139         args.push_back( make_pair( sizeof(cl_int) , (void *)&contour->total ));
 140         args.push_back( make_pair( sizeof(cl_mem) , (void *)&reader_oclmat.data ));
 141         args.push_back( make_pair( sizeof(cl_mem) , (void *)&dst_a.data ));
 142         cl_int dst_step = (cl_int)dst_a.step;
 143         args.push_back( make_pair( sizeof(cl_int) , (void *)&dst_step ));
 144
 145         openCLExecuteKernel(dst_a.clCxt, &moments, "icvContourMoments", globalThreads, localThreads, args, -1, -1);
 146
 147         cv::Mat dst(dst_a);
 148         a00 = a10 = a01 = a20 = a11 = a02 = a30 = a21 = a12 = a03 = 0.0;
 149         if (!cv::ocl::Context::getContext()->supportsFeature(Context::CL_DOUBLE))
 150         {
 151             for (int i = 0; i < contour->total; ++i)
 152             {
 153                 a00 += dst.at<cl_long>(0, i);
 154                 a10 += dst.at<cl_long>(1, i);
 155                 a01 += dst.at<cl_long>(2, i);
 156                 a20 += dst.at<cl_long>(3, i);
 157                 a11 += dst.at<cl_long>(4, i);
 158                 a02 += dst.at<cl_long>(5, i);
 159                 a30 += dst.at<cl_long>(6, i);
 160                 a21 += dst.at<cl_long>(7, i);
 161                 a12 += dst.at<cl_long>(8, i);
 162                 a03 += dst.at<cl_long>(9, i);
 163             }
 164         }
 165         else
 166         {
 167             a00 = cv::sum(dst.row(0))[0];
 168             a10 = cv::sum(dst.row(1))[0];
 169             a01 = cv::sum(dst.row(2))[0];
 170             a20 = cv::sum(dst.row(3))[0];
 171             a11 = cv::sum(dst.row(4))[0];
 172             a02 = cv::sum(dst.row(5))[0];
 173             a30 = cv::sum(dst.row(6))[0];
 174             a21 = cv::sum(dst.row(7))[0];
 175             a12 = cv::sum(dst.row(8))[0];
 176             a03 = cv::sum(dst.row(9))[0];
 177         }
 178
 179         double db1_2, db1_6, db1_12, db1_24, db1_20, db1_60;
 180         if( fabs(a00) > FLT_EPSILON )
 181         {
 182             if( a00 > 0 )
 183             {
 184                 db1_2 = 0.5;
 185                 db1_6 = 0.16666666666666666666666666666667;
 186                 db1_12 = 0.083333333333333333333333333333333;
 187                 db1_24 = 0.041666666666666666666666666666667;
 188                 db1_20 = 0.05;
 189                 db1_60 = 0.016666666666666666666666666666667;
 190             }
 191             else
 192             {
 193                 db1_2 = -0.5;
 194                 db1_6 = -0.16666666666666666666666666666667;
 195                 db1_12 = -0.083333333333333333333333333333333;
 196                 db1_24 = -0.041666666666666666666666666666667;
 197                 db1_20 = -0.05;
 198                 db1_60 = -0.016666666666666666666666666666667;
 199             }
 200
 201             // spatial moments
 202             mom->m00 = a00 * db1_2;
 203             mom->m10 = a10 * db1_6;
 204             mom->m01 = a01 * db1_6;
 205             mom->m20 = a20 * db1_12;
 206             mom->m11 = a11 * db1_24;
 207             mom->m02 = a02 * db1_12;
 208             mom->m30 = a30 * db1_20;
 209             mom->m21 = a21 * db1_60;
 210             mom->m12 = a12 * db1_60;
 211             mom->m03 = a03 * db1_20;
 212
 213             icvCompleteMomentState( mom );
 214         }
 215     }
 216 }
 217
 218 static void ocl_cvMoments( const void* array, CvMoments* mom, int binary )
 219 {
 220     const int TILE_SIZE = 256;
 221     int type, depth, cn, coi = 0;
 222     CvMat stub, *mat = (CvMat*)array;
 223     CvContour contourHeader;
 224     CvSeq* contour = 0;
 225     CvSeqBlock block;
 226     if( CV_IS_SEQ( array ))
 227     {
 228         contour = (CvSeq*)array;
 229         if( !CV_IS_SEQ_POINT_SET( contour ))
 230             CV_Error( CV_StsBadArg, "The passed sequence is not a valid contour" );
 231     }
 232
 233     if( !moments )
 234         CV_Error( CV_StsNullPtr, "" );
 235
 236     memset( mom, 0, sizeof(*mom));
 237
 238     if( !contour )
 239     {
 240
 241         mat = cvGetMat( mat, &stub, &coi );
 242         type = CV_MAT_TYPE( mat->type );
 243
 244         if( type == CV_32SC2 || type == CV_32FC2 )
 245         {
 246             contour = cvPointSeqFromMat(
 247                           CV_SEQ_KIND_CURVE | CV_SEQ_FLAG_CLOSED,
 248                           mat, &contourHeader, &block );
 249         }
 250     }
 251     if( contour )
 252     {
 253         icvContourMoments( contour, mom );
 254         return;
 255     }
 256
 257     type = CV_MAT_TYPE( mat->type );
 258     depth = CV_MAT_DEPTH( type );
 259     cn = CV_MAT_CN( type );
 260
 261     cv::Size size = cvGetMatSize( mat );
 262     if( cn > 1 && coi == 0 )
 263         CV_Error( CV_StsBadArg, "Invalid image type" );
 264
 265     if( size.width <= 0 || size.height <= 0 )
 266         return;
 267
 268     cv::Mat src0(mat);
 269     cv::ocl::oclMat src(src0);
 270     cv::Size tileSize;
 271     int blockx,blocky;
 272     if(size.width%TILE_SIZE == 0)
 273         blockx = size.width/TILE_SIZE;
 274     else
 275         blockx = size.width/TILE_SIZE + 1;
 276     if(size.height%TILE_SIZE == 0)
 277         blocky = size.height/TILE_SIZE;
 278     else
 279         blocky = size.height/TILE_SIZE + 1;
 280     oclMat dst_m(blocky * 10, blockx, CV_64FC1);
 281     oclMat sum(1, 10, CV_64FC1);
 282     int tile_width  = std::min(size.width,TILE_SIZE);
 283     int tile_height = std::min(size.height,TILE_SIZE);
 284     size_t localThreads[3]  = { tile_height, 1, 1};
 285     size_t globalThreads[3] = { size.height, blockx, 1};
 286     vector<pair<size_t , const void *> > args,args_sum;
 287     args.push_back( make_pair( sizeof(cl_mem) , (void *)&src.data ));
 288     args.push_back( make_pair( sizeof(cl_int) , (void *)&src.rows ));
 289     args.push_back( make_pair( sizeof(cl_int) , (void *)&src.cols ));
 290     args.push_back( make_pair( sizeof(cl_int) , (void *)&src.step ));
 291     args.push_back( make_pair( sizeof(cl_mem) , (void *)&dst_m.data ));
 292     args.push_back( make_pair( sizeof(cl_int) , (void *)&dst_m.cols ));
 293     args.push_back( make_pair( sizeof(cl_int) , (void *)&dst_m.step ));
 294     args.push_back( make_pair( sizeof(cl_int) , (void *)&blocky ));
 295     args.push_back( make_pair( sizeof(cl_int) , (void *)&depth ));
 296     args.push_back( make_pair( sizeof(cl_int) , (void *)&cn ));
 297     args.push_back( make_pair( sizeof(cl_int) , (void *)&coi ));
 298     args.push_back( make_pair( sizeof(cl_int) , (void *)&binary ));
 299     args.push_back( make_pair( sizeof(cl_int) , (void *)&TILE_SIZE ));
 300     openCLExecuteKernel(Context::getContext(), &moments, "CvMoments", globalThreads, localThreads, args, -1, depth);
 301
 302     size_t localThreadss[3]  = { 128, 1, 1};
 303     size_t globalThreadss[3] = { 128, 1, 1};
 304     args_sum.push_back( make_pair( sizeof(cl_int) , (void *)&src.rows ));
 305     args_sum.push_back( make_pair( sizeof(cl_int) , (void *)&src.cols ));
 306     args_sum.push_back( make_pair( sizeof(cl_int) , (void *)&tile_height ));
 307     args_sum.push_back( make_pair( sizeof(cl_int) , (void *)&tile_width ));
 308     args_sum.push_back( make_pair( sizeof(cl_int) , (void *)&TILE_SIZE ));
 309     args_sum.push_back( make_pair( sizeof(cl_mem) , (void *)&sum.data ));
 310     args_sum.push_back( make_pair( sizeof(cl_mem) , (void *)&dst_m.data ));
 311     args_sum.push_back( make_pair( sizeof(cl_int) , (void *)&dst_m.step ));
 312     openCLExecuteKernel(Context::getContext(), &moments, "dst_sum", globalThreadss, localThreadss, args_sum, -1, -1);
 313
 314     Mat dstsum(sum);
 315     mom->m00 = dstsum.at<double>(0, 0);
 316     mom->m10 = dstsum.at<double>(0, 1);
 317     mom->m01 = dstsum.at<double>(0, 2);
 318     mom->m20 = dstsum.at<double>(0, 3);
 319     mom->m11 = dstsum.at<double>(0, 4);
 320     mom->m02 = dstsum.at<double>(0, 5);
 321     mom->m30 = dstsum.at<double>(0, 6);
 322     mom->m21 = dstsum.at<double>(0, 7);
 323     mom->m12 = dstsum.at<double>(0, 8);
 324     mom->m03 = dstsum.at<double>(0, 9);
 325
 326     icvCompleteMomentState( mom );
 327 }
 328
 329 Moments ocl_moments( InputArray _array, bool binaryImage )
 330 {
 331     CvMoments om;
 332     Mat arr = _array.getMat();
 333     CvMat c_array = arr;
 334     ocl_cvMoments(&c_array, &om, binaryImage);
 335     return om;
 336 }
 337
 338 }
 339
 340 }
 341