moments.cpp

00001 /*M///////////////////////////////////////////////////////////////////////////////////////
00002 //
00003 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
00004 //
00005 //  By downloading, copying, installing or using the software you agree to this license.
00006 //  If you do not agree to this license, do not download, install,
00007 //  copy or use the software.
00008 //
00009 //
00010 //                        Intel License Agreement
00011 //                For Open Source Computer Vision Library
00012 //
00013 // Copyright (C) 2000, Intel Corporation, all rights reserved.
00014 // Third party copyrights are property of their respective owners.
00015 //
00016 // Redistribution and use in source and binary forms, with or without modification,
00017 // are permitted provided that the following conditions are met:
00018 //
00019 //   * Redistribution's of source code must retain the above copyright notice,
00020 //     this list of conditions and the following disclaimer.
00021 //
00022 //   * Redistribution's in binary form must reproduce the above copyright notice,
00023 //     this list of conditions and the following disclaimer in the documentation
00024 //     and/or other materials provided with the distribution.
00025 //
00026 //   * The name of Intel Corporation may not be used to endorse or promote products
00027 //     derived from this software without specific prior written permission.
00028 //
00029 // This software is provided by the copyright holders and contributors "as is" and
00030 // any express or implied warranties, including, but not limited to, the implied
00031 // warranties of merchantability and fitness for a particular purpose are disclaimed.
00032 // In no event shall the Intel Corporation or contributors be liable for any direct,
00033 // indirect, incidental, special, exemplary, or consequential damages
00034 // (including, but not limited to, procurement of substitute goods or services;
00035 // loss of use, data, or profits; or business interruption) however caused
00036 // and on any theory of liability, whether in contract, strict liability,
00037 // or tort (including negligence or otherwise) arising in any way out of
00038 // the use of this software, even if advised of the possibility of such damage.
00039 //
00040 //M*/
00041 #include "precomp.hpp"
00042 #include "opencl_kernels_imgproc.hpp"
00043 
00044 namespace cv
00045 {
00046 
00047 // The function calculates center of gravity and the central second order moments
00048 static void completeMomentState( Moments* moments )
00049 {
00050     double cx = 0, cy = 0;
00051     double mu20, mu11, mu02;
00052     double inv_m00 = 0.0;
00053     assert( moments != 0 );
00054 
00055     if( fabs(moments->m00) > DBL_EPSILON )
00056     {
00057         inv_m00 = 1. / moments->m00;
00058         cx = moments->m10 * inv_m00;
00059         cy = moments->m01 * inv_m00;
00060     }
00061 
00062     // mu20 = m20 - m10*cx
00063     mu20 = moments->m20 - moments->m10 * cx;
00064     // mu11 = m11 - m10*cy
00065     mu11 = moments->m11 - moments->m10 * cy;
00066     // mu02 = m02 - m01*cy
00067     mu02 = moments->m02 - moments->m01 * cy;
00068 
00069     moments->mu20 = mu20;
00070     moments->mu11 = mu11;
00071     moments->mu02 = mu02;
00072 
00073     // mu30 = m30 - cx*(3*mu20 + cx*m10)
00074     moments->mu30 = moments->m30 - cx * (3 * mu20 + cx * moments->m10);
00075     mu11 += mu11;
00076     // mu21 = m21 - cx*(2*mu11 + cx*m01) - cy*mu20
00077     moments->mu21 = moments->m21 - cx * (mu11 + cx * moments->m01) - cy * mu20;
00078     // mu12 = m12 - cy*(2*mu11 + cy*m10) - cx*mu02
00079     moments->mu12 = moments->m12 - cy * (mu11 + cy * moments->m10) - cx * mu02;
00080     // mu03 = m03 - cy*(3*mu02 + cy*m01)
00081     moments->mu03 = moments->m03 - cy * (3 * mu02 + cy * moments->m01);
00082 
00083 
00084     double inv_sqrt_m00 = std::sqrt(std::abs(inv_m00));
00085     double s2 = inv_m00*inv_m00, s3 = s2*inv_sqrt_m00;
00086 
00087     moments->nu20 = moments->mu20*s2; moments->nu11 = moments->mu11*s2; moments->nu02 = moments->mu02*s2;
00088     moments->nu30 = moments->mu30*s3; moments->nu21 = moments->mu21*s3; moments->nu12 = moments->mu12*s3; moments->nu03 = moments->mu03*s3;
00089 
00090 }
00091 
00092 
00093 static Moments contourMoments( const Mat& contour )
00094 {
00095     Moments m;
00096     int lpt = contour.checkVector(2);
00097     int is_float = contour.depth() == CV_32F;
00098     const Point* ptsi = contour.ptr<Point>();
00099     const Point2f* ptsf = contour.ptr<Point2f>();
00100 
00101     CV_Assert( contour.depth() == CV_32S || contour.depth() == CV_32F );
00102 
00103     if( lpt == 0 )
00104         return m;
00105 
00106     double a00 = 0, a10 = 0, a01 = 0, a20 = 0, a11 = 0, a02 = 0, a30 = 0, a21 = 0, a12 = 0, a03 = 0;
00107     double xi, yi, xi2, yi2, xi_1, yi_1, xi_12, yi_12, dxy, xii_1, yii_1;
00108 
00109     if( !is_float )
00110     {
00111         xi_1 = ptsi[lpt-1].x;
00112         yi_1 = ptsi[lpt-1].y;
00113     }
00114     else
00115     {
00116         xi_1 = ptsf[lpt-1].x;
00117         yi_1 = ptsf[lpt-1].y;
00118     }
00119 
00120     xi_12 = xi_1 * xi_1;
00121     yi_12 = yi_1 * yi_1;
00122 
00123     for( int i = 0; i < lpt; i++ )
00124     {
00125         if( !is_float )
00126         {
00127             xi = ptsi[i].x;
00128             yi = ptsi[i].y;
00129         }
00130         else
00131         {
00132             xi = ptsf[i].x;
00133             yi = ptsf[i].y;
00134         }
00135 
00136         xi2 = xi * xi;
00137         yi2 = yi * yi;
00138         dxy = xi_1 * yi - xi * yi_1;
00139         xii_1 = xi_1 + xi;
00140         yii_1 = yi_1 + yi;
00141 
00142         a00 += dxy;
00143         a10 += dxy * xii_1;
00144         a01 += dxy * yii_1;
00145         a20 += dxy * (xi_1 * xii_1 + xi2);
00146         a11 += dxy * (xi_1 * (yii_1 + yi_1) + xi * (yii_1 + yi));
00147         a02 += dxy * (yi_1 * yii_1 + yi2);
00148         a30 += dxy * xii_1 * (xi_12 + xi2);
00149         a03 += dxy * yii_1 * (yi_12 + yi2);
00150         a21 += dxy * (xi_12 * (3 * yi_1 + yi) + 2 * xi * xi_1 * yii_1 +
00151                    xi2 * (yi_1 + 3 * yi));
00152         a12 += dxy * (yi_12 * (3 * xi_1 + xi) + 2 * yi * yi_1 * xii_1 +
00153                    yi2 * (xi_1 + 3 * xi));
00154         xi_1 = xi;
00155         yi_1 = yi;
00156         xi_12 = xi2;
00157         yi_12 = yi2;
00158     }
00159 
00160     if( fabs(a00) > FLT_EPSILON )
00161     {
00162         double db1_2, db1_6, db1_12, db1_24, db1_20, db1_60;
00163 
00164         if( a00 > 0 )
00165         {
00166             db1_2 = 0.5;
00167             db1_6 = 0.16666666666666666666666666666667;
00168             db1_12 = 0.083333333333333333333333333333333;
00169             db1_24 = 0.041666666666666666666666666666667;
00170             db1_20 = 0.05;
00171             db1_60 = 0.016666666666666666666666666666667;
00172         }
00173         else
00174         {
00175             db1_2 = -0.5;
00176             db1_6 = -0.16666666666666666666666666666667;
00177             db1_12 = -0.083333333333333333333333333333333;
00178             db1_24 = -0.041666666666666666666666666666667;
00179             db1_20 = -0.05;
00180             db1_60 = -0.016666666666666666666666666666667;
00181         }
00182 
00183         // spatial moments
00184         m.m00 = a00 * db1_2;
00185         m.m10 = a10 * db1_6;
00186         m.m01 = a01 * db1_6;
00187         m.m20 = a20 * db1_12;
00188         m.m11 = a11 * db1_24;
00189         m.m02 = a02 * db1_12;
00190         m.m30 = a30 * db1_20;
00191         m.m21 = a21 * db1_60;
00192         m.m12 = a12 * db1_60;
00193         m.m03 = a03 * db1_20;
00194 
00195         completeMomentState( &m );
00196     }
00197     return m;
00198 }
00199 
00200 
00201 /****************************************************************************************\
00202 *                                Spatial Raster Moments                                  *
00203 \****************************************************************************************/
00204 
00205 template<typename T, typename WT, typename MT>
00206 struct MomentsInTile_SIMD
00207 {
00208     int operator() (const T *, int, WT &, WT &, WT &, MT &)
00209     {
00210         return 0;
00211     }
00212 };
00213 
00214 #if CV_SSE2
00215 
00216 template <>
00217 struct MomentsInTile_SIMD<uchar, int, int>
00218 {
00219     MomentsInTile_SIMD()
00220     {
00221         useSIMD = checkHardwareSupport(CV_CPU_SSE2);
00222     }
00223 
00224     int operator() (const uchar * ptr, int len, int & x0, int & x1, int & x2, int & x3)
00225     {
00226         int x = 0;
00227 
00228         if( useSIMD )
00229         {
00230             __m128i qx_init = _mm_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7);
00231             __m128i dx = _mm_set1_epi16(8);
00232             __m128i z = _mm_setzero_si128(), qx0 = z, qx1 = z, qx2 = z, qx3 = z, qx = qx_init;
00233 
00234             for( ; x <= len - 8; x += 8 )
00235             {
00236                 __m128i p = _mm_unpacklo_epi8(_mm_loadl_epi64((const __m128i*)(ptr + x)), z);
00237                 __m128i sx = _mm_mullo_epi16(qx, qx);
00238 
00239                 qx0 = _mm_add_epi32(qx0, _mm_sad_epu8(p, z));
00240                 qx1 = _mm_add_epi32(qx1, _mm_madd_epi16(p, qx));
00241                 qx2 = _mm_add_epi32(qx2, _mm_madd_epi16(p, sx));
00242                 qx3 = _mm_add_epi32(qx3, _mm_madd_epi16( _mm_mullo_epi16(p, qx), sx));
00243 
00244                 qx = _mm_add_epi16(qx, dx);
00245             }
00246 
00247             _mm_store_si128((__m128i*)buf, qx0);
00248             x0 = buf[0] + buf[1] + buf[2] + buf[3];
00249             _mm_store_si128((__m128i*)buf, qx1);
00250             x1 = buf[0] + buf[1] + buf[2] + buf[3];
00251             _mm_store_si128((__m128i*)buf, qx2);
00252             x2 = buf[0] + buf[1] + buf[2] + buf[3];
00253             _mm_store_si128((__m128i*)buf, qx3);
00254             x3 = buf[0] + buf[1] + buf[2] + buf[3];
00255         }
00256 
00257         return x;
00258     }
00259 
00260     int CV_DECL_ALIGNED(16) buf[4];
00261     bool useSIMD;
00262 };
00263 
00264 #elif CV_NEON
00265 
00266 template <>
00267 struct MomentsInTile_SIMD<uchar, int, int>
00268 {
00269     MomentsInTile_SIMD()
00270     {
00271         ushort CV_DECL_ALIGNED(8) init[4] = { 0, 1, 2, 3 };
00272         qx_init = vld1_u16(init);
00273         v_step = vdup_n_u16(4);
00274     }
00275 
00276     int operator() (const uchar * ptr, int len, int & x0, int & x1, int & x2, int & x3)
00277     {
00278         int x = 0;
00279 
00280         uint32x4_t v_z = vdupq_n_u32(0), v_x0 = v_z, v_x1 = v_z,
00281             v_x2 = v_z, v_x3 = v_z;
00282         uint16x4_t qx = qx_init;
00283 
00284         for( ; x <= len - 8; x += 8 )
00285         {
00286             uint16x8_t v_src = vmovl_u8(vld1_u8(ptr + x));
00287 
00288             // first part
00289             uint32x4_t v_qx = vmovl_u16(qx);
00290             uint16x4_t v_p = vget_low_u16(v_src);
00291             uint32x4_t v_px = vmull_u16(qx, v_p);
00292 
00293             v_x0 = vaddw_u16(v_x0, v_p);
00294             v_x1 = vaddq_u32(v_x1, v_px);
00295             v_px = vmulq_u32(v_px, v_qx);
00296             v_x2 = vaddq_u32(v_x2, v_px);
00297             v_x3 = vaddq_u32(v_x3, vmulq_u32(v_px, v_qx));
00298             qx = vadd_u16(qx, v_step);
00299 
00300             // second part
00301             v_qx = vmovl_u16(qx);
00302             v_p = vget_high_u16(v_src);
00303             v_px = vmull_u16(qx, v_p);
00304 
00305             v_x0 = vaddw_u16(v_x0, v_p);
00306             v_x1 = vaddq_u32(v_x1, v_px);
00307             v_px = vmulq_u32(v_px, v_qx);
00308             v_x2 = vaddq_u32(v_x2, v_px);
00309             v_x3 = vaddq_u32(v_x3, vmulq_u32(v_px, v_qx));
00310 
00311             qx = vadd_u16(qx, v_step);
00312         }
00313 
00314         vst1q_u32(buf, v_x0);
00315         x0 = buf[0] + buf[1] + buf[2] + buf[3];
00316         vst1q_u32(buf, v_x1);
00317         x1 = buf[0] + buf[1] + buf[2] + buf[3];
00318         vst1q_u32(buf, v_x2);
00319         x2 = buf[0] + buf[1] + buf[2] + buf[3];
00320         vst1q_u32(buf, v_x3);
00321         x3 = buf[0] + buf[1] + buf[2] + buf[3];
00322 
00323         return x;
00324     }
00325 
00326     uint CV_DECL_ALIGNED(16) buf[4];
00327     uint16x4_t qx_init, v_step;
00328 };
00329 
00330 #endif
00331 
00332 #if CV_SSE4_1
00333 
00334 template <>
00335 struct MomentsInTile_SIMD<ushort, int, int64>
00336 {
00337     MomentsInTile_SIMD()
00338     {
00339         useSIMD = checkHardwareSupport(CV_CPU_SSE4_1);
00340     }
00341 
00342     int operator() (const ushort * ptr, int len, int & x0, int & x1, int & x2, int64 & x3)
00343     {
00344         int x = 0;
00345 
00346         if (useSIMD)
00347         {
00348             __m128i vx_init0 = _mm_setr_epi32(0, 1, 2, 3), vx_init1 = _mm_setr_epi32(4, 5, 6, 7),
00349                 v_delta = _mm_set1_epi32(8), v_zero = _mm_setzero_si128(), v_x0 = v_zero,
00350                 v_x1 = v_zero, v_x2 = v_zero, v_x3 = v_zero, v_ix0 = vx_init0, v_ix1 = vx_init1;
00351 
00352             for( ; x <= len - 8; x += 8 )
00353             {
00354                 __m128i v_src = _mm_loadu_si128((const __m128i *)(ptr + x));
00355                 __m128i v_src0 = _mm_unpacklo_epi16(v_src, v_zero), v_src1 = _mm_unpackhi_epi16(v_src, v_zero);
00356 
00357                 v_x0 = _mm_add_epi32(v_x0, _mm_add_epi32(v_src0, v_src1));
00358                 __m128i v_x1_0 = _mm_mullo_epi32(v_src0, v_ix0), v_x1_1 = _mm_mullo_epi32(v_src1, v_ix1);
00359                 v_x1 = _mm_add_epi32(v_x1, _mm_add_epi32(v_x1_0, v_x1_1));
00360 
00361                 __m128i v_2ix0 = _mm_mullo_epi32(v_ix0, v_ix0), v_2ix1 = _mm_mullo_epi32(v_ix1, v_ix1);
00362                 v_x2 = _mm_add_epi32(v_x2, _mm_add_epi32(_mm_mullo_epi32(v_2ix0, v_src0), _mm_mullo_epi32(v_2ix1, v_src1)));
00363 
00364                 __m128i t = _mm_add_epi32(_mm_mullo_epi32(v_2ix0, v_x1_0), _mm_mullo_epi32(v_2ix1, v_x1_1));
00365                 v_x3 = _mm_add_epi64(v_x3, _mm_add_epi64(_mm_unpacklo_epi32(t, v_zero), _mm_unpackhi_epi32(t, v_zero)));
00366 
00367                 v_ix0 = _mm_add_epi32(v_ix0, v_delta);
00368                 v_ix1 = _mm_add_epi32(v_ix1, v_delta);
00369             }
00370 
00371             _mm_store_si128((__m128i*)buf, v_x0);
00372             x0 = buf[0] + buf[1] + buf[2] + buf[3];
00373             _mm_store_si128((__m128i*)buf, v_x1);
00374             x1 = buf[0] + buf[1] + buf[2] + buf[3];
00375             _mm_store_si128((__m128i*)buf, v_x2);
00376             x2 = buf[0] + buf[1] + buf[2] + buf[3];
00377 
00378             _mm_store_si128((__m128i*)buf64, v_x3);
00379             x3 = buf64[0] + buf64[1];
00380         }
00381 
00382         return x;
00383     }
00384 
00385     int CV_DECL_ALIGNED(16) buf[4];
00386     int64 CV_DECL_ALIGNED(16) buf64[2];
00387     bool useSIMD;
00388 };
00389 
00390 #endif
00391 
00392 template<typename T, typename WT, typename MT>
00393 #if defined __GNUC__ && __GNUC__ == 4 && __GNUC_MINOR__ >= 5 && __GNUC_MINOR__ < 9
00394 // Workaround for http://gcc.gnu.org/bugzilla/show_bug.cgi?id=60196
00395 __attribute__((optimize("no-tree-vectorize")))
00396 #endif
00397 static void momentsInTile( const Mat& img, double* moments )
00398 {
00399     Size size = img.size();
00400     int x, y;
00401     MT mom[10] = {0,0,0,0,0,0,0,0,0,0};
00402     MomentsInTile_SIMD<T, WT, MT> vop;
00403 
00404     for( y = 0; y < size.height; y++ )
00405     {
00406         const T* ptr = img.ptr<T>(y);
00407         WT x0 = 0, x1 = 0, x2 = 0;
00408         MT x3 = 0;
00409         x = vop(ptr, size.width, x0, x1, x2, x3);
00410 
00411         for( ; x < size.width; x++ )
00412         {
00413             WT p = ptr[x];
00414             WT xp = x * p, xxp;
00415 
00416             x0 += p;
00417             x1 += xp;
00418             xxp = xp * x;
00419             x2 += xxp;
00420             x3 += xxp * x;
00421         }
00422 
00423         WT py = y * x0, sy = y*y;
00424 
00425         mom[9] += ((MT)py) * sy;  // m03
00426         mom[8] += ((MT)x1) * sy;  // m12
00427         mom[7] += ((MT)x2) * y;  // m21
00428         mom[6] += x3;             // m30
00429         mom[5] += x0 * sy;        // m02
00430         mom[4] += x1 * y;         // m11
00431         mom[3] += x2;             // m20
00432         mom[2] += py;             // m01
00433         mom[1] += x1;             // m10
00434         mom[0] += x0;             // m00
00435     }
00436 
00437     for( x = 0; x < 10; x++ )
00438         moments[x] = (double)mom[x];
00439 }
00440 
00441 typedef void (*MomentsInTileFunc)(const Mat& img, double* moments);
00442 
00443 Moments::Moments()
00444 {
00445     m00 = m10 = m01 = m20 = m11 = m02 = m30 = m21 = m12 = m03 =
00446     mu20 = mu11 = mu02 = mu30 = mu21 = mu12 = mu03 =
00447     nu20 = nu11 = nu02 = nu30 = nu21 = nu12 = nu03 = 0.;
00448 }
00449 
00450 Moments::Moments( double _m00, double _m10, double _m01, double _m20, double _m11,
00451                   double _m02, double _m30, double _m21, double _m12, double _m03 )
00452 {
00453     m00 = _m00; m10 = _m10; m01 = _m01;
00454     m20 = _m20; m11 = _m11; m02 = _m02;
00455     m30 = _m30; m21 = _m21; m12 = _m12; m03 = _m03;
00456 
00457     double cx = 0, cy = 0, inv_m00 = 0;
00458     if( std::abs(m00) > DBL_EPSILON )
00459     {
00460         inv_m00 = 1./m00;
00461         cx = m10*inv_m00; cy = m01*inv_m00;
00462     }
00463 
00464     mu20 = m20 - m10*cx;
00465     mu11 = m11 - m10*cy;
00466     mu02 = m02 - m01*cy;
00467 
00468     mu30 = m30 - cx*(3*mu20 + cx*m10);
00469     mu21 = m21 - cx*(2*mu11 + cx*m01) - cy*mu20;
00470     mu12 = m12 - cy*(2*mu11 + cy*m10) - cx*mu02;
00471     mu03 = m03 - cy*(3*mu02 + cy*m01);
00472 
00473     double inv_sqrt_m00 = std::sqrt(std::abs(inv_m00));
00474     double s2 = inv_m00*inv_m00, s3 = s2*inv_sqrt_m00;
00475 
00476     nu20 = mu20*s2; nu11 = mu11*s2; nu02 = mu02*s2;
00477     nu30 = mu30*s3; nu21 = mu21*s3; nu12 = mu12*s3; nu03 = mu03*s3;
00478 }
00479 
00480 #ifdef HAVE_OPENCL
00481 
00482 static bool ocl_moments( InputArray _src, Moments& m, bool binary)
00483 {
00484     const int TILE_SIZE = 32;
00485     const int K = 10;
00486 
00487     ocl::Kernel k = ocl::Kernel("moments", ocl::imgproc::moments_oclsrc,
00488         format("-D TILE_SIZE=%d%s",
00489         TILE_SIZE,
00490         binary ? " -D OP_MOMENTS_BINARY" : ""));
00491 
00492     if( k.empty() )
00493         return false;
00494 
00495     UMat  src = _src.getUMat();
00496     Size sz = src.size();
00497     int xtiles = (sz.width + TILE_SIZE-1)/TILE_SIZE;
00498     int ytiles = (sz.height + TILE_SIZE-1)/TILE_SIZE;
00499     int ntiles = xtiles*ytiles;
00500     UMat  umbuf(1, ntiles*K, CV_32S);
00501 
00502     size_t globalsize[] = {(size_t)xtiles, (size_t)sz.height}, localsize[] = {1, TILE_SIZE};
00503     bool ok = k.args(ocl::KernelArg::ReadOnly(src),
00504                      ocl::KernelArg::PtrWriteOnly(umbuf),
00505                      xtiles).run(2, globalsize, localsize, true);
00506     if(!ok)
00507         return false;
00508     Mat mbuf = umbuf.getMat(ACCESS_READ);
00509     for( int i = 0; i < ntiles; i++ )
00510     {
00511         double x = (i % xtiles)*TILE_SIZE, y = (i / xtiles)*TILE_SIZE;
00512         const int* mom = mbuf.ptr<int>() + i*K;
00513         double xm = x * mom[0], ym = y * mom[0];
00514 
00515         // accumulate moments computed in each tile
00516 
00517         // + m00 ( = m00' )
00518         m.m00 += mom[0];
00519 
00520         // + m10 ( = m10' + x*m00' )
00521         m.m10 += mom[1] + xm;
00522 
00523         // + m01 ( = m01' + y*m00' )
00524         m.m01 += mom[2] + ym;
00525 
00526         // + m20 ( = m20' + 2*x*m10' + x*x*m00' )
00527         m.m20 += mom[3] + x * (mom[1] * 2 + xm);
00528 
00529         // + m11 ( = m11' + x*m01' + y*m10' + x*y*m00' )
00530         m.m11 += mom[4] + x * (mom[2] + ym) + y * mom[1];
00531 
00532         // + m02 ( = m02' + 2*y*m01' + y*y*m00' )
00533         m.m02 += mom[5] + y * (mom[2] * 2 + ym);
00534 
00535         // + m30 ( = m30' + 3*x*m20' + 3*x*x*m10' + x*x*x*m00' )
00536         m.m30 += mom[6] + x * (3. * mom[3] + x * (3. * mom[1] + xm));
00537 
00538         // + m21 ( = m21' + x*(2*m11' + 2*y*m10' + x*m01' + x*y*m00') + y*m20')
00539         m.m21 += mom[7] + x * (2 * (mom[4] + y * mom[1]) + x * (mom[2] + ym)) + y * mom[3];
00540 
00541         // + m12 ( = m12' + y*(2*m11' + 2*x*m01' + y*m10' + x*y*m00') + x*m02')
00542         m.m12 += mom[8] + y * (2 * (mom[4] + x * mom[2]) + y * (mom[1] + xm)) + x * mom[5];
00543 
00544         // + m03 ( = m03' + 3*y*m02' + 3*y*y*m01' + y*y*y*m00' )
00545         m.m03 += mom[9] + y * (3. * mom[5] + y * (3. * mom[2] + ym));
00546     }
00547 
00548     return true;
00549 }
00550 
00551 #endif
00552 
00553 }
00554 
00555 
00556 cv::Moments cv::moments( InputArray _src, bool binary )
00557 {
00558     const int TILE_SIZE = 32;
00559     MomentsInTileFunc func = 0;
00560     uchar nzbuf[TILE_SIZE*TILE_SIZE];
00561     Moments m;
00562     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
00563     Size size = _src.size();
00564 
00565     if( size.width <= 0 || size.height <= 0 )
00566         return m;
00567 
00568 #ifdef HAVE_OPENCL
00569     if( !(ocl::useOpenCL() && type == CV_8UC1  &&
00570         _src.isUMat() && ocl_moments(_src, m, binary)) )
00571 #endif
00572     {
00573         Mat mat = _src.getMat();
00574         if( mat.checkVector(2) >= 0 && (depth == CV_32F || depth == CV_32S))
00575             return contourMoments(mat);
00576 
00577         if( cn > 1 )
00578             CV_Error( CV_StsBadArg, "Invalid image type (must be single-channel)" );
00579 
00580 #if IPP_VERSION_X100 >= 810 && IPP_DISABLE_BLOCK
00581         CV_IPP_CHECK()
00582         {
00583             if (!binary)
00584             {
00585                 IppiSize roi = { mat.cols, mat.rows };
00586                 IppiMomentState_64f * moment = NULL;
00587                 // ippiMomentInitAlloc_64f, ippiMomentFree_64f are deprecated in 8.1, but there are not another way
00588                 // to initialize IppiMomentState_64f. When GetStateSize and Init functions will appear we have to
00589                 // change our code.
00590                 CV_SUPPRESS_DEPRECATED_START
00591                 if (ippiMomentInitAlloc_64f(&moment, ippAlgHintAccurate) >= 0)
00592                 {
00593                     typedef IppStatus (CV_STDCALL * ippiMoments)(const void * pSrc, int srcStep, IppiSize roiSize, IppiMomentState_64f* pCtx);
00594                     ippiMoments ippFunc =
00595                         type == CV_8UC1 ? (ippiMoments)ippiMoments64f_8u_C1R :
00596                         type == CV_16UC1 ? (ippiMoments)ippiMoments64f_16u_C1R :
00597                         type == CV_32FC1? (ippiMoments)ippiMoments64f_32f_C1R : 0;
00598 
00599                     if (ippFunc)
00600                     {
00601                         if (ippFunc(mat.data, (int)mat.step, roi, moment) >= 0)
00602                         {
00603                             IppiPoint point = { 0, 0 };
00604                             ippiGetSpatialMoment_64f(moment, 0, 0, 0, point, &m.m00);
00605                             ippiGetSpatialMoment_64f(moment, 1, 0, 0, point, &m.m10);
00606                             ippiGetSpatialMoment_64f(moment, 0, 1, 0, point, &m.m01);
00607 
00608                             ippiGetSpatialMoment_64f(moment, 2, 0, 0, point, &m.m20);
00609                             ippiGetSpatialMoment_64f(moment, 1, 1, 0, point, &m.m11);
00610                             ippiGetSpatialMoment_64f(moment, 0, 2, 0, point, &m.m02);
00611 
00612                             ippiGetSpatialMoment_64f(moment, 3, 0, 0, point, &m.m30);
00613                             ippiGetSpatialMoment_64f(moment, 2, 1, 0, point, &m.m21);
00614                             ippiGetSpatialMoment_64f(moment, 1, 2, 0, point, &m.m12);
00615                             ippiGetSpatialMoment_64f(moment, 0, 3, 0, point, &m.m03);
00616                             ippiGetCentralMoment_64f(moment, 2, 0, 0, &m.mu20);
00617                             ippiGetCentralMoment_64f(moment, 1, 1, 0, &m.mu11);
00618                             ippiGetCentralMoment_64f(moment, 0, 2, 0, &m.mu02);
00619                             ippiGetCentralMoment_64f(moment, 3, 0, 0, &m.mu30);
00620                             ippiGetCentralMoment_64f(moment, 2, 1, 0, &m.mu21);
00621                             ippiGetCentralMoment_64f(moment, 1, 2, 0, &m.mu12);
00622                             ippiGetCentralMoment_64f(moment, 0, 3, 0, &m.mu03);
00623                             ippiGetNormalizedCentralMoment_64f(moment, 2, 0, 0, &m.nu20);
00624                             ippiGetNormalizedCentralMoment_64f(moment, 1, 1, 0, &m.nu11);
00625                             ippiGetNormalizedCentralMoment_64f(moment, 0, 2, 0, &m.nu02);
00626                             ippiGetNormalizedCentralMoment_64f(moment, 3, 0, 0, &m.nu30);
00627                             ippiGetNormalizedCentralMoment_64f(moment, 2, 1, 0, &m.nu21);
00628                             ippiGetNormalizedCentralMoment_64f(moment, 1, 2, 0, &m.nu12);
00629                             ippiGetNormalizedCentralMoment_64f(moment, 0, 3, 0, &m.nu03);
00630 
00631                             ippiMomentFree_64f(moment);
00632                             CV_IMPL_ADD(CV_IMPL_IPP);
00633                             return m;
00634                         }
00635                         setIppErrorStatus();
00636                     }
00637                     ippiMomentFree_64f(moment);
00638                 }
00639                 else
00640                     setIppErrorStatus();
00641                 CV_SUPPRESS_DEPRECATED_END
00642             }
00643         }
00644 #endif
00645 
00646         if( binary || depth == CV_8U )
00647             func = momentsInTile<uchar, int, int>;
00648         else if( depth == CV_16U )
00649             func = momentsInTile<ushort, int, int64>;
00650         else if( depth == CV_16S )
00651             func = momentsInTile<short, int, int64>;
00652         else if( depth == CV_32F )
00653             func = momentsInTile<float, double, double>;
00654         else if( depth == CV_64F )
00655             func = momentsInTile<double, double, double>;
00656         else
00657             CV_Error( CV_StsUnsupportedFormat, "" );
00658 
00659         Mat src0(mat);
00660 
00661         for( int y = 0; y < size.height; y += TILE_SIZE )
00662         {
00663             Size tileSize;
00664             tileSize.height = std::min(TILE_SIZE, size.height - y);
00665 
00666             for( int x = 0; x < size.width; x += TILE_SIZE )
00667             {
00668                 tileSize.width = std::min(TILE_SIZE, size.width - x);
00669                 Mat src(src0, cv::Rect(x, y, tileSize.width, tileSize.height));
00670 
00671                 if( binary )
00672                 {
00673                     cv::Mat tmp(tileSize, CV_8U, nzbuf);
00674                     cv::compare( src, 0, tmp, CV_CMP_NE );
00675                     src = tmp;
00676                 }
00677 
00678                 double mom[10];
00679                 func( src, mom );
00680 
00681                 if(binary)
00682                 {
00683                     double s = 1./255;
00684                     for( int k = 0; k < 10; k++ )
00685                         mom[k] *= s;
00686                 }
00687 
00688                 double xm = x * mom[0], ym = y * mom[0];
00689 
00690                 // accumulate moments computed in each tile
00691 
00692                 // + m00 ( = m00' )
00693                 m.m00 += mom[0];
00694 
00695                 // + m10 ( = m10' + x*m00' )
00696                 m.m10 += mom[1] + xm;
00697 
00698                 // + m01 ( = m01' + y*m00' )
00699                 m.m01 += mom[2] + ym;
00700 
00701                 // + m20 ( = m20' + 2*x*m10' + x*x*m00' )
00702                 m.m20 += mom[3] + x * (mom[1] * 2 + xm);
00703 
00704                 // + m11 ( = m11' + x*m01' + y*m10' + x*y*m00' )
00705                 m.m11 += mom[4] + x * (mom[2] + ym) + y * mom[1];
00706 
00707                 // + m02 ( = m02' + 2*y*m01' + y*y*m00' )
00708                 m.m02 += mom[5] + y * (mom[2] * 2 + ym);
00709 
00710                 // + m30 ( = m30' + 3*x*m20' + 3*x*x*m10' + x*x*x*m00' )
00711                 m.m30 += mom[6] + x * (3. * mom[3] + x * (3. * mom[1] + xm));
00712 
00713                 // + m21 ( = m21' + x*(2*m11' + 2*y*m10' + x*m01' + x*y*m00') + y*m20')
00714                 m.m21 += mom[7] + x * (2 * (mom[4] + y * mom[1]) + x * (mom[2] + ym)) + y * mom[3];
00715 
00716                 // + m12 ( = m12' + y*(2*m11' + 2*x*m01' + y*m10' + x*y*m00') + x*m02')
00717                 m.m12 += mom[8] + y * (2 * (mom[4] + x * mom[2]) + y * (mom[1] + xm)) + x * mom[5];
00718 
00719                 // + m03 ( = m03' + 3*y*m02' + 3*y*y*m01' + y*y*y*m00' )
00720                 m.m03 += mom[9] + y * (3. * mom[5] + y * (3. * mom[2] + ym));
00721             }
00722         }
00723     }
00724 
00725     completeMomentState( &m );
00726     return m;
00727 }
00728 
00729 
00730 void cv::HuMoments( const Moments& m, double hu[7] )
00731 {
00732     double t0 = m.nu30 + m.nu12;
00733     double t1 = m.nu21 + m.nu03;
00734 
00735     double q0 = t0 * t0, q1 = t1 * t1;
00736 
00737     double n4 = 4 * m.nu11;
00738     double s = m.nu20 + m.nu02;
00739     double d = m.nu20 - m.nu02;
00740 
00741     hu[0] = s;
00742     hu[1] = d * d + n4 * m.nu11;
00743     hu[3] = q0 + q1;
00744     hu[5] = d * (q0 - q1) + n4 * t0 * t1;
00745 
00746     t0 *= q0 - 3 * q1;
00747     t1 *= 3 * q0 - q1;
00748 
00749     q0 = m.nu30 - 3 * m.nu12;
00750     q1 = 3 * m.nu21 - m.nu03;
00751 
00752     hu[2] = q0 * q0 + q1 * q1;
00753     hu[4] = q0 * t0 + q1 * t1;
00754     hu[6] = q1 * t0 - q0 * t1;
00755 }
00756 
00757 void cv::HuMoments( const Moments& m, OutputArray _hu )
00758 {
00759     _hu.create(7, 1, CV_64F);
00760     Mat hu = _hu.getMat();
00761     CV_Assert( hu.isContinuous() );
00762     HuMoments(m, hu.ptr<double>());
00763 }
00764 
00765 
00766 CV_IMPL void cvMoments( const CvArr* arr, CvMoments* moments, int binary )
00767 {
00768     const IplImage* img = (const IplImage*)arr;
00769     cv::Mat src;
00770     if( CV_IS_IMAGE(arr) && img->roi && img->roi->coi > 0 )
00771         cv::extractImageCOI(arr, src, img->roi->coi-1);
00772     else
00773         src = cv::cvarrToMat(arr);
00774     cv::Moments m = cv::moments(src, binary != 0);
00775     CV_Assert( moments != 0 );
00776     *moments = m;
00777 }
00778 
00779 
00780 CV_IMPL double cvGetSpatialMoment( CvMoments * moments, int x_order, int y_order )
00781 {
00782     int order = x_order + y_order;
00783 
00784     if( !moments )
00785         CV_Error( CV_StsNullPtr, "" );
00786     if( (x_order | y_order) < 0 || order > 3 )
00787         CV_Error( CV_StsOutOfRange, "" );
00788 
00789     return (&(moments->m00))[order + (order >> 1) + (order > 2) * 2 + y_order];
00790 }
00791 
00792 
00793 CV_IMPL double cvGetCentralMoment( CvMoments * moments, int x_order, int y_order )
00794 {
00795     int order = x_order + y_order;
00796 
00797     if( !moments )
00798         CV_Error( CV_StsNullPtr, "" );
00799     if( (x_order | y_order) < 0 || order > 3 )
00800         CV_Error( CV_StsOutOfRange, "" );
00801 
00802     return order >= 2 ? (&(moments->m00))[4 + order * 3 + y_order] :
00803     order == 0 ? moments->m00 : 0;
00804 }
00805 
00806 
00807 CV_IMPL double cvGetNormalizedCentralMoment( CvMoments * moments, int x_order, int y_order )
00808 {
00809     int order = x_order + y_order;
00810 
00811     double mu = cvGetCentralMoment( moments, x_order, y_order );
00812     double m00s = moments->inv_sqrt_m00;
00813 
00814     while( --order >= 0 )
00815         mu *= m00s;
00816     return mu * m00s * m00s;
00817 }
00818 
00819 
00820 CV_IMPL void cvGetHuMoments( CvMoments * mState, CvHuMoments * HuState )
00821 {
00822     if( !mState || !HuState )
00823         CV_Error( CV_StsNullPtr, "" );
00824 
00825     double m00s = mState->inv_sqrt_m00, m00 = m00s * m00s, s2 = m00 * m00, s3 = s2 * m00s;
00826 
00827     double nu20 = mState->mu20 * s2,
00828     nu11 = mState->mu11 * s2,
00829     nu02 = mState->mu02 * s2,
00830     nu30 = mState->mu30 * s3,
00831     nu21 = mState->mu21 * s3, nu12 = mState->mu12 * s3, nu03 = mState->mu03 * s3;
00832 
00833     double t0 = nu30 + nu12;
00834     double t1 = nu21 + nu03;
00835 
00836     double q0 = t0 * t0, q1 = t1 * t1;
00837 
00838     double n4 = 4 * nu11;
00839     double s = nu20 + nu02;
00840     double d = nu20 - nu02;
00841 
00842     HuState->hu1 = s;
00843     HuState->hu2 = d * d + n4 * nu11;
00844     HuState->hu4 = q0 + q1;
00845     HuState->hu6 = d * (q0 - q1) + n4 * t0 * t1;
00846 
00847     t0 *= q0 - 3 * q1;
00848     t1 *= 3 * q0 - q1;
00849 
00850     q0 = nu30 - 3 * nu12;
00851     q1 = 3 * nu21 - nu03;
00852 
00853     HuState->hu3 = q0 * q0 + q1 * q1;
00854     HuState->hu5 = q0 * t0 + q1 * t1;
00855     HuState->hu7 = q1 * t0 - q0 * t1;
00856 }
00857 
00858 
00859 /* End of file. */
00860