modules/imgproc/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 //                        Intel License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
  14 // Third party copyrights are property of their respective owners.
  15 //
  16 // Redistribution and use in source and binary forms, with or without modification,
  17 // are permitted provided that the following conditions are met:
  18 //
  19 //   * Redistribution's of source code must retain the above copyright notice,
  20 //     this list of conditions and the following disclaimer.
  21 //
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,
  23 //     this list of conditions and the following disclaimer in the documentation
  24 //     and/or other materials provided with the distribution.
  25 //
  26 //   * The name of Intel Corporation may not be used to endorse or promote products
  27 //     derived from this software without specific prior written permission.
  28 //
  29 // This software is provided by the copyright holders and contributors "as is" and
  30 // any express or implied warranties, including, but not limited to, the implied
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,
  33 // indirect, incidental, special, exemplary, or consequential damages
  34 // (including, but not limited to, procurement of substitute goods or services;
  35 // loss of use, data, or profits; or business interruption) however caused
  36 // and on any theory of liability, whether in contract, strict liability,
  37 // or tort (including negligence or otherwise) arising in any way out of
  38 // the use of this software, even if advised of the possibility of such damage.
  39 //
  40 //M*/
  41 #include "precomp.hpp"
  42
  43 // The function calculates center of gravity and the central second order moments
  44 static void icvCompleteMomentState( CvMoments* moments )
  45 {
  46     double cx = 0, cy = 0;
  47     double mu20, mu11, mu02;
  48
  49     assert( moments != 0 );
  50     moments->inv_sqrt_m00 = 0;
  51
  52     if( fabs(moments->m00) > DBL_EPSILON )
  53     {
  54         double inv_m00 = 1. / moments->m00;
  55         cx = moments->m10 * inv_m00;
  56         cy = moments->m01 * inv_m00;
  57         moments->inv_sqrt_m00 = std::sqrt( fabs(inv_m00) );
  58     }
  59
  60     // mu20 = m20 - m10*cx
  61     mu20 = moments->m20 - moments->m10 * cx;
  62     // mu11 = m11 - m10*cy
  63     mu11 = moments->m11 - moments->m10 * cy;
  64     // mu02 = m02 - m01*cy
  65     mu02 = moments->m02 - moments->m01 * cy;
  66
  67     moments->mu20 = mu20;
  68     moments->mu11 = mu11;
  69     moments->mu02 = mu02;
  70
  71     // mu30 = m30 - cx*(3*mu20 + cx*m10)
  72     moments->mu30 = moments->m30 - cx * (3 * mu20 + cx * moments->m10);
  73     mu11 += mu11;
  74     // mu21 = m21 - cx*(2*mu11 + cx*m01) - cy*mu20
  75     moments->mu21 = moments->m21 - cx * (mu11 + cx * moments->m01) - cy * mu20;
  76     // mu12 = m12 - cy*(2*mu11 + cy*m10) - cx*mu02
  77     moments->mu12 = moments->m12 - cy * (mu11 + cy * moments->m10) - cx * mu02;
  78     // mu03 = m03 - cy*(3*mu02 + cy*m01)
  79     moments->mu03 = moments->m03 - cy * (3 * mu02 + cy * moments->m01);
  80 }
  81
  82
  83 static void icvContourMoments( CvSeq* contour, CvMoments* moments )
  84 {
  85     int is_float = CV_SEQ_ELTYPE(contour) == CV_32FC2;
  86
  87     if( contour->total )
  88     {
  89         CvSeqReader reader;
  90         double a00, a10, a01, a20, a11, a02, a30, a21, a12, a03;
  91         double xi, yi, xi2, yi2, xi_1, yi_1, xi_12, yi_12, dxy, xii_1, yii_1;
  92         int lpt = contour->total;
  93
  94         a00 = a10 = a01 = a20 = a11 = a02 = a30 = a21 = a12 = a03 = 0;
  95
  96         cvStartReadSeq( contour, &reader, 0 );
  97
  98         if( !is_float )
  99         {
 100             xi_1 = ((CvPoint*)(reader.ptr))->x;
 101             yi_1 = ((CvPoint*)(reader.ptr))->y;
 102         }
 103         else
 104         {
 105             xi_1 = ((CvPoint2D32f*)(reader.ptr))->x;
 106             yi_1 = ((CvPoint2D32f*)(reader.ptr))->y;
 107         }
 108         CV_NEXT_SEQ_ELEM( contour->elem_size, reader );
 109
 110         xi_12 = xi_1 * xi_1;
 111         yi_12 = yi_1 * yi_1;
 112
 113         while( lpt-- > 0 )
 114         {
 115             if( !is_float )
 116             {
 117                 xi = ((CvPoint*)(reader.ptr))->x;
 118                 yi = ((CvPoint*)(reader.ptr))->y;
 119             }
 120             else
 121             {
 122                 xi = ((CvPoint2D32f*)(reader.ptr))->x;
 123                 yi = ((CvPoint2D32f*)(reader.ptr))->y;
 124             }
 125             CV_NEXT_SEQ_ELEM( contour->elem_size, reader );
 126
 127             xi2 = xi * xi;
 128             yi2 = yi * yi;
 129             dxy = xi_1 * yi - xi * yi_1;
 130             xii_1 = xi_1 + xi;
 131             yii_1 = yi_1 + yi;
 132
 133             a00 += dxy;
 134             a10 += dxy * xii_1;
 135             a01 += dxy * yii_1;
 136             a20 += dxy * (xi_1 * xii_1 + xi2);
 137             a11 += dxy * (xi_1 * (yii_1 + yi_1) + xi * (yii_1 + yi));
 138             a02 += dxy * (yi_1 * yii_1 + yi2);
 139             a30 += dxy * xii_1 * (xi_12 + xi2);
 140             a03 += dxy * yii_1 * (yi_12 + yi2);
 141             a21 +=
 142                 dxy * (xi_12 * (3 * yi_1 + yi) + 2 * xi * xi_1 * yii_1 +
 143                        xi2 * (yi_1 + 3 * yi));
 144             a12 +=
 145                 dxy * (yi_12 * (3 * xi_1 + xi) + 2 * yi * yi_1 * xii_1 +
 146                        yi2 * (xi_1 + 3 * xi));
 147
 148             xi_1 = xi;
 149             yi_1 = yi;
 150             xi_12 = xi2;
 151             yi_12 = yi2;
 152         }
 153
 154         double db1_2, db1_6, db1_12, db1_24, db1_20, db1_60;
 155
 156         if( fabs(a00) > FLT_EPSILON )
 157         {
 158             if( a00 > 0 )
 159             {
 160                 db1_2 = 0.5;
 161                 db1_6 = 0.16666666666666666666666666666667;
 162                 db1_12 = 0.083333333333333333333333333333333;
 163                 db1_24 = 0.041666666666666666666666666666667;
 164                 db1_20 = 0.05;
 165                 db1_60 = 0.016666666666666666666666666666667;
 166             }
 167             else
 168             {
 169                 db1_2 = -0.5;
 170                 db1_6 = -0.16666666666666666666666666666667;
 171                 db1_12 = -0.083333333333333333333333333333333;
 172                 db1_24 = -0.041666666666666666666666666666667;
 173                 db1_20 = -0.05;
 174                 db1_60 = -0.016666666666666666666666666666667;
 175             }
 176
 177             // spatial moments
 178             moments->m00 = a00 * db1_2;
 179             moments->m10 = a10 * db1_6;
 180             moments->m01 = a01 * db1_6;
 181             moments->m20 = a20 * db1_12;
 182             moments->m11 = a11 * db1_24;
 183             moments->m02 = a02 * db1_12;
 184             moments->m30 = a30 * db1_20;
 185             moments->m21 = a21 * db1_60;
 186             moments->m12 = a12 * db1_60;
 187             moments->m03 = a03 * db1_20;
 188
 189             icvCompleteMomentState( moments );
 190         }
 191     }
 192 }
 193
 194
 195 /****************************************************************************************\
 196 *                                Spatial Raster Moments                                  *
 197 \****************************************************************************************/
 198
 199 template<typename T, typename WT, typename MT>
 200 #if defined __GNUC__ && __GNUC__ == 4 && __GNUC_MINOR__ >= 5 && __GNUC_MINOR__ < 9
 201 // Workaround for http://gcc.gnu.org/bugzilla/show_bug.cgi?id=60196
 202 __attribute__((optimize("no-tree-vectorize")))
 203 #endif
 204 static void momentsInTile( const cv::Mat& img, double* moments )
 205 {
 206     cv::Size size = img.size();
 207     int x, y;
 208     MT mom[10] = {0,0,0,0,0,0,0,0,0,0};
 209
 210     for( y = 0; y < size.height; y++ )
 211     {
 212         const T* ptr = (const T*)(img.data + y*img.step);
 213         WT x0 = 0, x1 = 0, x2 = 0;
 214         MT x3 = 0;
 215
 216         for( x = 0; x < size.width; x++ )
 217         {
 218             WT p = ptr[x];
 219             WT xp = x * p, xxp;
 220
 221             x0 += p;
 222             x1 += xp;
 223             xxp = xp * x;
 224             x2 += xxp;
 225             x3 += xxp * x;
 226         }
 227
 228         WT py = y * x0, sy = y*y;
 229
 230         mom[9] += ((MT)py) * sy;  // m03
 231         mom[8] += ((MT)x1) * sy;  // m12
 232         mom[7] += ((MT)x2) * y;  // m21
 233         mom[6] += x3;             // m30
 234         mom[5] += x0 * sy;        // m02
 235         mom[4] += x1 * y;         // m11
 236         mom[3] += x2;             // m20
 237         mom[2] += py;             // m01
 238         mom[1] += x1;             // m10
 239         mom[0] += x0;             // m00
 240     }
 241
 242     for( x = 0; x < 10; x++ )
 243         moments[x] = (double)mom[x];
 244 }
 245
 246
 247 #if CV_SSE2
 248
 249 template<> void momentsInTile<uchar, int, int>( const cv::Mat& img, double* moments )
 250 {
 251     typedef uchar T;
 252     typedef int WT;
 253     typedef int MT;
 254     cv::Size size = img.size();
 255     int y;
 256     MT mom[10] = {0,0,0,0,0,0,0,0,0,0};
 257     bool useSIMD = cv::checkHardwareSupport(CV_CPU_SSE2);
 258
 259     for( y = 0; y < size.height; y++ )
 260     {
 261         const T* ptr = img.ptr<T>(y);
 262         int x0 = 0, x1 = 0, x2 = 0, x3 = 0, x = 0;
 263
 264         if( useSIMD )
 265         {
 266             __m128i qx_init = _mm_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7);
 267             __m128i dx = _mm_set1_epi16(8);
 268             __m128i z = _mm_setzero_si128(), qx0 = z, qx1 = z, qx2 = z, qx3 = z, qx = qx_init;
 269
 270             for( ; x <= size.width - 8; x += 8 )
 271             {
 272                 __m128i p = _mm_unpacklo_epi8(_mm_loadl_epi64((const __m128i*)(ptr + x)), z);
 273                 qx0 = _mm_add_epi32(qx0, _mm_sad_epu8(p, z));
 274                 __m128i px = _mm_mullo_epi16(p, qx);
 275                 __m128i sx = _mm_mullo_epi16(qx, qx);
 276                 qx1 = _mm_add_epi32(qx1, _mm_madd_epi16(p, qx));
 277                 qx2 = _mm_add_epi32(qx2, _mm_madd_epi16(p, sx));
 278                 qx3 = _mm_add_epi32(qx3, _mm_madd_epi16(px, sx));
 279
 280                 qx = _mm_add_epi16(qx, dx);
 281             }
 282             int CV_DECL_ALIGNED(16) buf[4];
 283             _mm_store_si128((__m128i*)buf, qx0);
 284             x0 = buf[0] + buf[1] + buf[2] + buf[3];
 285             _mm_store_si128((__m128i*)buf, qx1);
 286             x1 = buf[0] + buf[1] + buf[2] + buf[3];
 287             _mm_store_si128((__m128i*)buf, qx2);
 288             x2 = buf[0] + buf[1] + buf[2] + buf[3];
 289             _mm_store_si128((__m128i*)buf, qx3);
 290             x3 = buf[0] + buf[1] + buf[2] + buf[3];
 291         }
 292
 293         for( ; x < size.width; x++ )
 294         {
 295             WT p = ptr[x];
 296             WT xp = x * p, xxp;
 297
 298             x0 += p;
 299             x1 += xp;
 300             xxp = xp * x;
 301             x2 += xxp;
 302             x3 += xxp * x;
 303         }
 304
 305         WT py = y * x0, sy = y*y;
 306
 307         mom[9] += ((MT)py) * sy;  // m03
 308         mom[8] += ((MT)x1) * sy;  // m12
 309         mom[7] += ((MT)x2) * y;  // m21
 310         mom[6] += x3;             // m30
 311         mom[5] += x0 * sy;        // m02
 312         mom[4] += x1 * y;         // m11
 313         mom[3] += x2;             // m20
 314         mom[2] += py;             // m01
 315         mom[1] += x1;             // m10
 316         mom[0] += x0;             // m00
 317     }
 318
 319     for(int x = 0; x < 10; x++ )
 320         moments[x] = (double)mom[x];
 321 }
 322
 323 #endif
 324
 325 typedef void (*CvMomentsInTileFunc)(const cv::Mat& img, double* moments);
 326
 327 CV_IMPL void cvMoments( const void* array, CvMoments* moments, int binary )
 328 {
 329     const int TILE_SIZE = 32;
 330     int type, depth, cn, coi = 0;
 331     CvMat stub, *mat = (CvMat*)array;
 332     CvMomentsInTileFunc func = 0;
 333     CvContour contourHeader;
 334     CvSeq* contour = 0;
 335     CvSeqBlock block;
 336     double buf[TILE_SIZE*TILE_SIZE];
 337     uchar nzbuf[TILE_SIZE*TILE_SIZE];
 338
 339     if( CV_IS_SEQ( array ))
 340     {
 341         contour = (CvSeq*)array;
 342         if( !CV_IS_SEQ_POINT_SET( contour ))
 343             CV_Error( CV_StsBadArg, "The passed sequence is not a valid contour" );
 344     }
 345
 346     if( !moments )
 347         CV_Error( CV_StsNullPtr, "" );
 348
 349     memset( moments, 0, sizeof(*moments));
 350
 351     if( !contour )
 352     {
 353         mat = cvGetMat( mat, &stub, &coi );
 354         type = CV_MAT_TYPE( mat->type );
 355
 356         if( type == CV_32SC2 || type == CV_32FC2 )
 357         {
 358             contour = cvPointSeqFromMat(
 359                 CV_SEQ_KIND_CURVE | CV_SEQ_FLAG_CLOSED,
 360                 mat, &contourHeader, &block );
 361         }
 362     }
 363
 364     if( contour )
 365     {
 366         icvContourMoments( contour, moments );
 367         return;
 368     }
 369
 370     type = CV_MAT_TYPE( mat->type );
 371     depth = CV_MAT_DEPTH( type );
 372     cn = CV_MAT_CN( type );
 373
 374     cv::Size size = cvGetMatSize( mat );
 375
 376     if( cn > 1 && coi == 0 )
 377         CV_Error( CV_StsBadArg, "Invalid image type" );
 378
 379     if( size.width <= 0 || size.height <= 0 )
 380         return;
 381
 382     if( binary || depth == CV_8U )
 383         func = momentsInTile<uchar, int, int>;
 384     else if( depth == CV_16U )
 385         func = momentsInTile<ushort, int, int64>;
 386     else if( depth == CV_16S )
 387         func = momentsInTile<short, int, int64>;
 388     else if( depth == CV_32F )
 389         func = momentsInTile<float, double, double>;
 390     else if( depth == CV_64F )
 391         func = momentsInTile<double, double, double>;
 392     else
 393         CV_Error( CV_StsUnsupportedFormat, "" );
 394
 395     cv::Mat src0(mat);
 396
 397     for( int y = 0; y < size.height; y += TILE_SIZE )
 398     {
 399         cv::Size tileSize;
 400         tileSize.height = std::min(TILE_SIZE, size.height - y);
 401
 402         for( int x = 0; x < size.width; x += TILE_SIZE )
 403         {
 404             tileSize.width = std::min(TILE_SIZE, size.width - x);
 405             cv::Mat src(src0, cv::Rect(x, y, tileSize.width, tileSize.height));
 406
 407             if( coi > 0 )
 408             {
 409                 cv::Mat tmp(tileSize, depth, buf);
 410                 int pairs[] = {coi-1, 0};
 411                 cv::mixChannels(&src, 1, &tmp, 1, pairs, 1);
 412                 src = tmp;
 413             }
 414             if( binary )
 415             {
 416                 cv::Mat tmp(tileSize, CV_8U, nzbuf);
 417                 cv::compare( src, 0, tmp, CV_CMP_NE );
 418                 src = tmp;
 419             }
 420
 421             double mom[10];
 422             func( src, mom );
 423
 424             if(binary)
 425             {
 426                 double s = 1./255;
 427                 for( int k = 0; k < 10; k++ )
 428                     mom[k] *= s;
 429             }
 430
 431             double xm = x * mom[0], ym = y * mom[0];
 432
 433             // accumulate moments computed in each tile
 434
 435             // + m00 ( = m00' )
 436             moments->m00 += mom[0];
 437
 438             // + m10 ( = m10' + x*m00' )
 439             moments->m10 += mom[1] + xm;
 440
 441             // + m01 ( = m01' + y*m00' )
 442             moments->m01 += mom[2] + ym;
 443
 444             // + m20 ( = m20' + 2*x*m10' + x*x*m00' )
 445             moments->m20 += mom[3] + x * (mom[1] * 2 + xm);
 446
 447             // + m11 ( = m11' + x*m01' + y*m10' + x*y*m00' )
 448             moments->m11 += mom[4] + x * (mom[2] + ym) + y * mom[1];
 449
 450             // + m02 ( = m02' + 2*y*m01' + y*y*m00' )
 451             moments->m02 += mom[5] + y * (mom[2] * 2 + ym);
 452
 453             // + m30 ( = m30' + 3*x*m20' + 3*x*x*m10' + x*x*x*m00' )
 454             moments->m30 += mom[6] + x * (3. * mom[3] + x * (3. * mom[1] + xm));
 455
 456             // + m21 ( = m21' + x*(2*m11' + 2*y*m10' + x*m01' + x*y*m00') + y*m20')
 457             moments->m21 += mom[7] + x * (2 * (mom[4] + y * mom[1]) + x * (mom[2] + ym)) + y * mom[3];
 458
 459             // + m12 ( = m12' + y*(2*m11' + 2*x*m01' + y*m10' + x*y*m00') + x*m02')
 460             moments->m12 += mom[8] + y * (2 * (mom[4] + x * mom[2]) + y * (mom[1] + xm)) + x * mom[5];
 461
 462             // + m03 ( = m03' + 3*y*m02' + 3*y*y*m01' + y*y*y*m00' )
 463             moments->m03 += mom[9] + y * (3. * mom[5] + y * (3. * mom[2] + ym));
 464         }
 465     }
 466
 467     icvCompleteMomentState( moments );
 468 }
 469
 470
 471 CV_IMPL void cvGetHuMoments( CvMoments * mState, CvHuMoments * HuState )
 472 {
 473     if( !mState || !HuState )
 474         CV_Error( CV_StsNullPtr, "" );
 475
 476     double m00s = mState->inv_sqrt_m00, m00 = m00s * m00s, s2 = m00 * m00, s3 = s2 * m00s;
 477
 478     double nu20 = mState->mu20 * s2,
 479         nu11 = mState->mu11 * s2,
 480         nu02 = mState->mu02 * s2,
 481         nu30 = mState->mu30 * s3,
 482         nu21 = mState->mu21 * s3, nu12 = mState->mu12 * s3, nu03 = mState->mu03 * s3;
 483
 484     double t0 = nu30 + nu12;
 485     double t1 = nu21 + nu03;
 486
 487     double q0 = t0 * t0, q1 = t1 * t1;
 488
 489     double n4 = 4 * nu11;
 490     double s = nu20 + nu02;
 491     double d = nu20 - nu02;
 492
 493     HuState->hu1 = s;
 494     HuState->hu2 = d * d + n4 * nu11;
 495     HuState->hu4 = q0 + q1;
 496     HuState->hu6 = d * (q0 - q1) + n4 * t0 * t1;
 497
 498     t0 *= q0 - 3 * q1;
 499     t1 *= 3 * q0 - q1;
 500
 501     q0 = nu30 - 3 * nu12;
 502     q1 = 3 * nu21 - nu03;
 503
 504     HuState->hu3 = q0 * q0 + q1 * q1;
 505     HuState->hu5 = q0 * t0 + q1 * t1;
 506     HuState->hu7 = q1 * t0 - q0 * t1;
 507 }
 508
 509
 510 CV_IMPL double cvGetSpatialMoment( CvMoments * moments, int x_order, int y_order )
 511 {
 512     int order = x_order + y_order;
 513
 514     if( !moments )
 515         CV_Error( CV_StsNullPtr, "" );
 516     if( (x_order | y_order) < 0 || order > 3 )
 517         CV_Error( CV_StsOutOfRange, "" );
 518
 519     return (&(moments->m00))[order + (order >> 1) + (order > 2) * 2 + y_order];
 520 }
 521
 522
 523 CV_IMPL double cvGetCentralMoment( CvMoments * moments, int x_order, int y_order )
 524 {
 525     int order = x_order + y_order;
 526
 527     if( !moments )
 528         CV_Error( CV_StsNullPtr, "" );
 529     if( (x_order | y_order) < 0 || order > 3 )
 530         CV_Error( CV_StsOutOfRange, "" );
 531
 532     return order >= 2 ? (&(moments->m00))[4 + order * 3 + y_order] :
 533            order == 0 ? moments->m00 : 0;
 534 }
 535
 536
 537 CV_IMPL double cvGetNormalizedCentralMoment( CvMoments * moments, int x_order, int y_order )
 538 {
 539     int order = x_order + y_order;
 540
 541     double mu = cvGetCentralMoment( moments, x_order, y_order );
 542     double m00s = moments->inv_sqrt_m00;
 543
 544     while( --order >= 0 )
 545         mu *= m00s;
 546     return mu * m00s * m00s;
 547 }
 548
 549
 550 namespace cv
 551 {
 552
 553 Moments::Moments()
 554 {
 555     m00 = m10 = m01 = m20 = m11 = m02 = m30 = m21 = m12 = m03 =
 556     mu20 = mu11 = mu02 = mu30 = mu21 = mu12 = mu03 =
 557     nu20 = nu11 = nu02 = nu30 = nu21 = nu12 = nu03 = 0.;
 558 }
 559
 560 Moments::Moments( double _m00, double _m10, double _m01, double _m20, double _m11,
 561                   double _m02, double _m30, double _m21, double _m12, double _m03 )
 562 {
 563     m00 = _m00; m10 = _m10; m01 = _m01;
 564     m20 = _m20; m11 = _m11; m02 = _m02;
 565     m30 = _m30; m21 = _m21; m12 = _m12; m03 = _m03;
 566
 567     double cx = 0, cy = 0, inv_m00 = 0;
 568     if( std::abs(m00) > DBL_EPSILON )
 569     {
 570         inv_m00 = 1./m00;
 571         cx = m10*inv_m00; cy = m01*inv_m00;
 572     }
 573
 574     mu20 = m20 - m10*cx;
 575     mu11 = m11 - m10*cy;
 576     mu02 = m02 - m01*cy;
 577
 578     mu30 = m30 - cx*(3*mu20 + cx*m10);
 579     mu21 = m21 - cx*(2*mu11 + cx*m01) - cy*mu20;
 580     mu12 = m12 - cy*(2*mu11 + cy*m10) - cx*mu02;
 581     mu03 = m03 - cy*(3*mu02 + cy*m01);
 582
 583     double inv_sqrt_m00 = std::sqrt(std::abs(inv_m00));
 584     double s2 = inv_m00*inv_m00, s3 = s2*inv_sqrt_m00;
 585
 586     nu20 = mu20*s2; nu11 = mu11*s2; nu02 = mu02*s2;
 587     nu30 = mu30*s3; nu21 = mu21*s3; nu12 = mu12*s3; nu03 = mu03*s3;
 588 }
 589
 590 Moments::Moments( const CvMoments& m )
 591 {
 592     *this = Moments(m.m00, m.m10, m.m01, m.m20, m.m11, m.m02, m.m30, m.m21, m.m12, m.m03);
 593 }
 594
 595 Moments::operator CvMoments() const
 596 {
 597     CvMoments m;
 598     m.m00 = m00; m.m10 = m10; m.m01 = m01;
 599     m.m20 = m20; m.m11 = m11; m.m02 = m02;
 600     m.m30 = m30; m.m21 = m21; m.m12 = m12; m.m03 = m03;
 601     m.mu20 = mu20; m.mu11 = mu11; m.mu02 = mu02;
 602     m.mu30 = mu30; m.mu21 = mu21; m.mu12 = mu12; m.mu03 = mu03;
 603     double am00 = std::abs(m00);
 604     m.inv_sqrt_m00 = am00 > DBL_EPSILON ? 1./std::sqrt(am00) : 0;
 605
 606     return m;
 607 }
 608
 609 }
 610
 611 cv::Moments cv::moments( InputArray _array, bool binaryImage )
 612 {
 613     CvMoments om;
 614     Mat arr = _array.getMat();
 615     CvMat c_array = arr;
 616     cvMoments(&c_array, &om, binaryImage);
 617     return om;
 618 }
 619
 620 void cv::HuMoments( const Moments& m, double hu[7] )
 621 {
 622     double t0 = m.nu30 + m.nu12;
 623     double t1 = m.nu21 + m.nu03;
 624
 625     double q0 = t0 * t0, q1 = t1 * t1;
 626
 627     double n4 = 4 * m.nu11;
 628     double s = m.nu20 + m.nu02;
 629     double d = m.nu20 - m.nu02;
 630
 631     hu[0] = s;
 632     hu[1] = d * d + n4 * m.nu11;
 633     hu[3] = q0 + q1;
 634     hu[5] = d * (q0 - q1) + n4 * t0 * t1;
 635
 636     t0 *= q0 - 3 * q1;
 637     t1 *= 3 * q0 - q1;
 638
 639     q0 = m.nu30 - 3 * m.nu12;
 640     q1 = 3 * m.nu21 - m.nu03;
 641
 642     hu[2] = q0 * q0 + q1 * q1;
 643     hu[4] = q0 * t0 + q1 * t1;
 644     hu[6] = q1 * t0 - q0 * t1;
 645 }
 646
 647 void cv::HuMoments( const Moments& m, OutputArray _hu )
 648 {
 649     _hu.create(7, 1, CV_64F);
 650     Mat hu = _hu.getMat();
 651     CV_Assert( hu.isContinuous() );
 652     HuMoments(m, (double*)hu.data);
 653 }
 654
 655 /* End of file. */