git.maemo.org Git - opencv/blob - cv/src/cvmoments.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 "_cv.h"
  42
  43 /* The function calculates center of gravity and central second order moments */
  44 static void
  45 icvCompleteMomentState( CvMoments* moments )
  46 {
  47     double cx = 0, cy = 0;
  48     double mu20, mu11, mu02;
  49
  50     assert( moments != 0 );
  51     moments->inv_sqrt_m00 = 0;
  52
  53     if( fabs(moments->m00) > DBL_EPSILON )
  54     {
  55         double inv_m00 = 1. / moments->m00;
  56         cx = moments->m10 * inv_m00;
  57         cy = moments->m01 * inv_m00;
  58         moments->inv_sqrt_m00 = sqrt( fabs(inv_m00) );
  59     }
  60
  61     /* mu20 = m20 - m10*cx */
  62     mu20 = moments->m20 - moments->m10 * cx;
  63     /* mu11 = m11 - m10*cy */
  64     mu11 = moments->m11 - moments->m10 * cy;
  65     /* mu02 = m02 - m01*cy */
  66     mu02 = moments->m02 - moments->m01 * cy;
  67
  68     moments->mu20 = mu20;
  69     moments->mu11 = mu11;
  70     moments->mu02 = mu02;
  71
  72     /* mu30 = m30 - cx*(3*mu20 + cx*m10) */
  73     moments->mu30 = moments->m30 - cx * (3 * mu20 + cx * moments->m10);
  74     mu11 += mu11;
  75     /* mu21 = m21 - cx*(2*mu11 + cx*m01) - cy*mu20 */
  76     moments->mu21 = moments->m21 - cx * (mu11 + cx * moments->m01) - cy * mu20;
  77     /* mu12 = m12 - cy*(2*mu11 + cy*m10) - cx*mu02 */
  78     moments->mu12 = moments->m12 - cy * (mu11 + cy * moments->m10) - cx * mu02;
  79     /* mu03 = m03 - cy*(3*mu02 + cy*m01) */
  80     moments->mu03 = moments->m03 - cy * (3 * mu02 + cy * moments->m01);
  81 }
  82
  83
  84 static void
  85 icvContourMoments( CvSeq* contour, CvMoments* moments )
  86 {
  87     int is_float = CV_SEQ_ELTYPE(contour) == CV_32FC2;
  88
  89     if( contour->total )
  90     {
  91         CvSeqReader reader;
  92         double a00, a10, a01, a20, a11, a02, a30, a21, a12, a03;
  93         double xi, yi, xi2, yi2, xi_1, yi_1, xi_12, yi_12, dxy, xii_1, yii_1;
  94         int lpt = contour->total;
  95
  96         a00 = a10 = a01 = a20 = a11 = a02 = a30 = a21 = a12 = a03 = 0;
  97
  98         cvStartReadSeq( contour, &reader, 0 );
  99
 100         if( !is_float )
 101         {
 102             xi_1 = ((CvPoint*)(reader.ptr))->x;
 103             yi_1 = ((CvPoint*)(reader.ptr))->y;
 104         }
 105         else
 106         {
 107             xi_1 = ((CvPoint2D32f*)(reader.ptr))->x;
 108             yi_1 = ((CvPoint2D32f*)(reader.ptr))->y;
 109         }
 110         CV_NEXT_SEQ_ELEM( contour->elem_size, reader );
 111
 112         xi_12 = xi_1 * xi_1;
 113         yi_12 = yi_1 * yi_1;
 114
 115         while( lpt-- > 0 )
 116         {
 117             if( !is_float )
 118             {
 119                 xi = ((CvPoint*)(reader.ptr))->x;
 120                 yi = ((CvPoint*)(reader.ptr))->y;
 121             }
 122             else
 123             {
 124                 xi = ((CvPoint2D32f*)(reader.ptr))->x;
 125                 yi = ((CvPoint2D32f*)(reader.ptr))->y;
 126             }
 127             CV_NEXT_SEQ_ELEM( contour->elem_size, reader );
 128
 129             xi2 = xi * xi;
 130             yi2 = yi * yi;
 131             dxy = xi_1 * yi - xi * yi_1;
 132             xii_1 = xi_1 + xi;
 133             yii_1 = yi_1 + yi;
 134
 135             a00 += dxy;
 136             a10 += dxy * xii_1;
 137             a01 += dxy * yii_1;
 138             a20 += dxy * (xi_1 * xii_1 + xi2);
 139             a11 += dxy * (xi_1 * (yii_1 + yi_1) + xi * (yii_1 + yi));
 140             a02 += dxy * (yi_1 * yii_1 + yi2);
 141             a30 += dxy * xii_1 * (xi_12 + xi2);
 142             a03 += dxy * yii_1 * (yi_12 + yi2);
 143             a21 +=
 144                 dxy * (xi_12 * (3 * yi_1 + yi) + 2 * xi * xi_1 * yii_1 +
 145                        xi2 * (yi_1 + 3 * yi));
 146             a12 +=
 147                 dxy * (yi_12 * (3 * xi_1 + xi) + 2 * yi * yi_1 * xii_1 +
 148                        yi2 * (xi_1 + 3 * xi));
 149
 150             xi_1 = xi;
 151             yi_1 = yi;
 152             xi_12 = xi2;
 153             yi_12 = yi2;
 154         }
 155
 156         double db1_2, db1_6, db1_12, db1_24, db1_20, db1_60;
 157
 158         if( fabs(a00) > FLT_EPSILON )
 159         {
 160             if( a00 > 0 )
 161             {
 162                 db1_2 = 0.5;
 163                 db1_6 = 0.16666666666666666666666666666667;
 164                 db1_12 = 0.083333333333333333333333333333333;
 165                 db1_24 = 0.041666666666666666666666666666667;
 166                 db1_20 = 0.05;
 167                 db1_60 = 0.016666666666666666666666666666667;
 168             }
 169             else
 170             {
 171                 db1_2 = -0.5;
 172                 db1_6 = -0.16666666666666666666666666666667;
 173                 db1_12 = -0.083333333333333333333333333333333;
 174                 db1_24 = -0.041666666666666666666666666666667;
 175                 db1_20 = -0.05;
 176                 db1_60 = -0.016666666666666666666666666666667;
 177             }
 178
 179             /*  spatial moments    */
 180             moments->m00 = a00 * db1_2;
 181             moments->m10 = a10 * db1_6;
 182             moments->m01 = a01 * db1_6;
 183             moments->m20 = a20 * db1_12;
 184             moments->m11 = a11 * db1_24;
 185             moments->m02 = a02 * db1_12;
 186             moments->m30 = a30 * db1_20;
 187             moments->m21 = a21 * db1_60;
 188             moments->m12 = a12 * db1_60;
 189             moments->m03 = a03 * db1_20;
 190
 191             icvCompleteMomentState( moments );
 192         }
 193     }
 194 }
 195
 196
 197 /* summarizes moment values for all tiles */
 198 static void
 199 icvAccumulateMoments( double *tiles, CvSize size, CvSize tile_size, CvMoments * moments )
 200 {
 201     int x, y;
 202
 203     for( y = 0; y < size.height; y += tile_size.height )
 204     {
 205         for( x = 0; x < size.width; x += tile_size.width, tiles += 10 )
 206         {
 207             double dx = x, dy = y;
 208             double dxm = dx * tiles[0], dym = dy * tiles[0];
 209
 210             /* + m00 ( = m00' ) */
 211             moments->m00 += tiles[0];
 212
 213             /* + m10 ( = m10' + dx*m00' ) */
 214             moments->m10 += tiles[1] + dxm;
 215
 216             /* + m01 ( = m01' + dy*m00' ) */
 217             moments->m01 += tiles[2] + dym;
 218
 219             /* + m20 ( = m20' + 2*dx*m10' + dx*dx*m00' ) */
 220             moments->m20 += tiles[3] + dx * (tiles[1] * 2 + dxm);
 221
 222             /* + m11 ( = m11' + dx*m01' + dy*m10' + dx*dy*m00' ) */
 223             moments->m11 += tiles[4] + dx * (tiles[2] + dym) + dy * tiles[1];
 224
 225             /* + m02 ( = m02' + 2*dy*m01' + dy*dy*m00' ) */
 226             moments->m02 += tiles[5] + dy * (tiles[2] * 2 + dym);
 227
 228             /* + m30 ( = m30' + 3*dx*m20' + 3*dx*dx*m10' + dx*dx*dx*m00' ) */
 229             moments->m30 += tiles[6] + dx * (3. * tiles[3] + dx * (3. * tiles[1] + dxm));
 230
 231             /* + m21 (= m21' + dx*(2*m11' + 2*dy*m10' + dx*m01' + dx*dy*m00') + dy*m20') */
 232             moments->m21 += tiles[7] + dx * (2 * (tiles[4] + dy * tiles[1]) +
 233                                              dx * (tiles[2] + dym)) + dy * tiles[3];
 234
 235             /* + m12 (= m12' + dy*(2*m11' + 2*dx*m01' + dy*m10' + dx*dy*m00') + dx*m02') */
 236             moments->m12 += tiles[8] + dy * (2 * (tiles[4] + dx * tiles[2]) +
 237                                              dy * (tiles[1] + dxm)) + dx * tiles[5];
 238
 239             /* + m03 ( = m03' + 3*dy*m02' + 3*dy*dy*m01' + dy*dy*dy*m00' ) */
 240             moments->m03 += tiles[9] + dy * (3. * tiles[5] + dy * (3. * tiles[2] + dym));
 241         }
 242     }
 243
 244     icvCompleteMomentState( moments );
 245 }
 246
 247
 248 /****************************************************************************************\
 249 *                                   Spatial Moments                                      *
 250 \****************************************************************************************/
 251
 252 #define ICV_DEF_CALC_MOMENTS_IN_TILE( __op__, name, flavor, srctype, temptype, momtype ) \
 253 static CvStatus CV_STDCALL icv##name##_##flavor##_CnCR                                   \
 254 ( const srctype* img, int step, CvSize size, int cn, int coi, double *moments )          \
 255 {                                                                                        \
 256     int x, y, sx_init = (size.width & -4) * (size.width & -4), sy = 0;                   \
 257     momtype mom[10];                                                                     \
 258                                                                                          \
 259     assert( img && size.width && (size.width | size.height) >= 0 );                      \
 260     memset( mom, 0, 10 * sizeof( mom[0] ));                                              \
 261                                                                                          \
 262     if( coi )                                                                            \
 263         img += coi - 1;                                                                  \
 264     step /= sizeof(img[0]);                                                              \
 265                                                                                          \
 266     for( y = 0; y < size.height; sy += 2 * y + 1, y++, img += step )                     \
 267     {                                                                                    \
 268         temptype  x0 = 0;                                                                \
 269         temptype  x1 = 0;                                                                \
 270         temptype  x2 = 0;                                                                \
 271         momtype   x3 = 0;                                                                \
 272         int sx = sx_init;                                                                \
 273         const srctype* ptr = img;                                                        \
 274                                                                                          \
 275         for( x = 0; x < size.width - 3; x += 4, ptr += cn*4 )                            \
 276         {                                                                                \
 277             temptype p0 = __op__(ptr[0]), p1 = __op__(ptr[cn]),                          \
 278                      p2 = __op__(ptr[2*cn]), p3 = __op__(ptr[3*cn]);                     \
 279             temptype t = p1;                                                             \
 280             temptype a, b, c;                                                            \
 281                                                                                          \
 282             p0 += p1 + p2 + p3; /* p0 + p1 + p2 + p3 */                                  \
 283             p1 += 2 * p2 + 3 * p3;      /* p1 + p2*2 + p3*3 */                           \
 284             p2 = p1 + 2 * p2 + 6 * p3;  /* p1 + p2*4 + p3*9 */                           \
 285             p3 = 2 * p2 - t + 9 * p3;   /* p1 + p2*8 + p3*27 */                          \
 286                                                                                          \
 287             a = x * p0 + p1;    /* x*p0 + (x+1)*p1 + (x+2)*p2 + (x+3)*p3 */              \
 288             b = x * p1 + p2;    /* (x+1)*p1 + 2*(x+2)*p2 + 3*(x+3)*p3 */                 \
 289             c = x * p2 + p3;    /* (x+1)*p1 + 4*(x+2)*p2 + 9*(x+3)*p3 */                 \
 290                                                                                          \
 291             x0 += p0;                                                                    \
 292             x1 += a;                                                                     \
 293             a = a * x + b;      /*(x^2)*p0+((x+1)^2)*p1+((x+2)^2)*p2+((x+3)^2)*p3 */     \
 294             x2 += a;                                                                     \
 295             x3 += ((momtype)(a + b)) * x + c;  /*x3 += (x^3)*p0+((x+1)^3)*p1 +  */       \
 296                                                /*  ((x+2)^3)*p2+((x+3)^3)*p3   */        \
 297         }                                                                                \
 298                                                                                          \
 299         /* process the rest */                                                           \
 300         for( ; x < size.width; sx += 2 * x + 1, x++, ptr += cn )                         \
 301         {                                                                                \
 302             temptype p = __op__(ptr[0]);                                                 \
 303             temptype xp = x * p;                                                         \
 304                                                                                          \
 305             x0 += p;                                                                     \
 306             x1 += xp;                                                                    \
 307             x2 += sx * p;                                                                \
 308             x3 += ((momtype)sx) * xp;                                                    \
 309         }                                                                                \
 310                                                                                          \
 311         {                                                                                \
 312             temptype py = y * x0;                                                        \
 313                                                                                          \
 314             mom[9] += ((momtype)py) * sy;  /* m03 */                                     \
 315             mom[8] += ((momtype)x1) * sy;  /* m12 */                                     \
 316             mom[7] += ((momtype)x2) * y;   /* m21 */                                     \
 317             mom[6] += x3;                  /* m30 */                                     \
 318             mom[5] += x0 * sy;             /* m02 */                                     \
 319             mom[4] += x1 * y;              /* m11 */                                     \
 320             mom[3] += x2;                  /* m20 */                                     \
 321             mom[2] += py;                  /* m01 */                                     \
 322             mom[1] += x1;                  /* m10 */                                     \
 323             mom[0] += x0;                  /* m00 */                                     \
 324         }                                                                                \
 325     }                                                                                    \
 326                                                                                          \
 327     for( x = 0; x < 10; x++ )                                                            \
 328         moments[x] = (double)mom[x];                                                     \
 329                                                                                          \
 330     return CV_OK;                                                                        \
 331 }
 332
 333
 334 ICV_DEF_CALC_MOMENTS_IN_TILE( CV_NOP, MomentsInTile, 8u, uchar, int, int )
 335 ICV_DEF_CALC_MOMENTS_IN_TILE( CV_NOP, MomentsInTile, 16u, ushort, int, int64 )
 336 ICV_DEF_CALC_MOMENTS_IN_TILE( CV_NOP, MomentsInTile, 16s, short, int, int64 )
 337 ICV_DEF_CALC_MOMENTS_IN_TILE( CV_NOP, MomentsInTile, 32f, float, double, double )
 338 ICV_DEF_CALC_MOMENTS_IN_TILE( CV_NOP, MomentsInTile, 64f, double, double, double )
 339
 340 ICV_DEF_CALC_MOMENTS_IN_TILE( CV_NONZERO, MomentsInTileBin, 8u, uchar, int, int )
 341 ICV_DEF_CALC_MOMENTS_IN_TILE( CV_NONZERO, MomentsInTileBin, 16s, ushort, int, int )
 342 ICV_DEF_CALC_MOMENTS_IN_TILE( CV_NONZERO_FLT, MomentsInTileBin, 32f, int, int, int )
 343 ICV_DEF_CALC_MOMENTS_IN_TILE( CV_NONZERO_FLT, MomentsInTileBin, 64f, int64, double, double )
 344
 345 #define icvMomentsInTile_8s_CnCR  0
 346 #define icvMomentsInTile_32s_CnCR  0
 347 #define icvMomentsInTileBin_8s_CnCR   icvMomentsInTileBin_8u_CnCR
 348 #define icvMomentsInTileBin_16u_CnCR   icvMomentsInTileBin_16s_CnCR
 349 #define icvMomentsInTileBin_32s_CnCR  0
 350
 351 CV_DEF_INIT_FUNC_TAB_2D( MomentsInTile, CnCR )
 352 CV_DEF_INIT_FUNC_TAB_2D( MomentsInTileBin, CnCR )
 353
 354 ////////////////////////////////// IPP moment functions //////////////////////////////////
 355
 356 icvMoments_8u_C1R_t icvMoments_8u_C1R_p = 0;
 357 icvMoments_32f_C1R_t icvMoments_32f_C1R_p = 0;
 358 icvMomentInitAlloc_64f_t icvMomentInitAlloc_64f_p = 0;
 359 icvMomentFree_64f_t icvMomentFree_64f_p = 0;
 360 icvGetSpatialMoment_64f_t icvGetSpatialMoment_64f_p = 0;
 361
 362 typedef CvStatus (CV_STDCALL * CvMomentIPPFunc)
 363     ( const void* img, int step, CvSize size, void* momentstate );
 364
 365 CV_IMPL void
 366 cvMoments( const void* array, CvMoments* moments, int binary )
 367 {
 368     static CvFuncTable mom_tab;
 369     static CvFuncTable mombin_tab;
 370     static int inittab = 0;
 371     double* tiles = 0;
 372     void* ippmomentstate = 0;
 373
 374     CV_FUNCNAME("cvMoments");
 375
 376     __BEGIN__;
 377
 378     int type = 0, depth, cn, pix_size;
 379     int coi = 0;
 380     int x, y, k, tile_num = 1;
 381     CvSize size, tile_size = { 32, 32 };
 382     CvMat stub, *mat = (CvMat*)array;
 383     CvFunc2DnC_1A1P func = 0;
 384     CvMomentIPPFunc ipp_func = 0;
 385     CvContour contour_header;
 386     CvSeq* contour = 0;
 387     CvSeqBlock block;
 388
 389     if( CV_IS_SEQ( array ))
 390     {
 391         contour = (CvSeq*)array;
 392         if( !CV_IS_SEQ_POLYGON( contour ))
 393             CV_ERROR( CV_StsBadArg, "The passed sequence is not a valid contour" );
 394     }
 395
 396     if( !inittab )
 397     {
 398         icvInitMomentsInTileCnCRTable( &mom_tab );
 399         icvInitMomentsInTileBinCnCRTable( &mombin_tab );
 400         inittab = 1;
 401     }
 402
 403     if( !moments )
 404         CV_ERROR( CV_StsNullPtr, "" );
 405
 406     memset( moments, 0, sizeof(*moments));
 407
 408     if( !contour )
 409     {
 410         CV_CALL( mat = cvGetMat( mat, &stub, &coi ));
 411         type = CV_MAT_TYPE( mat->type );
 412
 413         if( type == CV_32SC2 || type == CV_32FC2 )
 414         {
 415             CV_CALL( contour = cvPointSeqFromMat(
 416                 CV_SEQ_KIND_CURVE | CV_SEQ_FLAG_CLOSED,
 417                 mat, &contour_header, &block ));
 418         }
 419     }
 420
 421     if( contour )
 422     {
 423         icvContourMoments( contour, moments );
 424         EXIT;
 425     }
 426
 427     type = CV_MAT_TYPE( mat->type );
 428     depth = CV_MAT_DEPTH( type );
 429     cn = CV_MAT_CN( type );
 430     pix_size = CV_ELEM_SIZE(type);
 431     size = cvGetMatSize( mat );
 432
 433     if( cn > 1 && coi == 0 )
 434         CV_ERROR( CV_StsBadArg, "Invalid image type" );
 435
 436     if( size.width <= 0 || size.height <= 0 )
 437     {
 438         EXIT;
 439     }
 440
 441     if( type == CV_8UC1 )
 442         ipp_func = (CvMomentIPPFunc)icvMoments_8u_C1R_p;
 443     else if( type == CV_32FC1 )
 444         ipp_func = (CvMomentIPPFunc)icvMoments_32f_C1R_p;
 445
 446     if( ipp_func && !binary )
 447     {
 448         int matstep = mat->step ? mat->step : CV_STUB_STEP;
 449         IPPI_CALL( icvMomentInitAlloc_64f_p( &ippmomentstate, cvAlgHintAccurate ));
 450         IPPI_CALL( ipp_func( mat->data.ptr, matstep, size, ippmomentstate ));
 451         icvGetSpatialMoment_64f_p( ippmomentstate, 0, 0, 0, cvPoint(0,0), &moments->m00 );
 452         icvGetSpatialMoment_64f_p( ippmomentstate, 1, 0, 0, cvPoint(0,0), &moments->m10 );
 453         icvGetSpatialMoment_64f_p( ippmomentstate, 0, 1, 0, cvPoint(0,0), &moments->m01 );
 454         icvGetSpatialMoment_64f_p( ippmomentstate, 2, 0, 0, cvPoint(0,0), &moments->m20 );
 455         icvGetSpatialMoment_64f_p( ippmomentstate, 1, 1, 0, cvPoint(0,0), &moments->m11 );
 456         icvGetSpatialMoment_64f_p( ippmomentstate, 0, 2, 0, cvPoint(0,0), &moments->m02 );
 457         icvGetSpatialMoment_64f_p( ippmomentstate, 3, 0, 0, cvPoint(0,0), &moments->m30 );
 458         icvGetSpatialMoment_64f_p( ippmomentstate, 2, 1, 0, cvPoint(0,0), &moments->m21 );
 459         icvGetSpatialMoment_64f_p( ippmomentstate, 1, 2, 0, cvPoint(0,0), &moments->m12 );
 460         icvGetSpatialMoment_64f_p( ippmomentstate, 0, 3, 0, cvPoint(0,0), &moments->m03 );
 461         icvCompleteMomentState( moments );
 462         EXIT;
 463     }
 464
 465     func = (CvFunc2DnC_1A1P)(!binary ? mom_tab.fn_2d[depth] : mombin_tab.fn_2d[depth]);
 466
 467     if( !func )
 468         CV_ERROR( CV_StsBadArg, cvUnsupportedFormat );
 469
 470     if( depth >= CV_32S && !binary )
 471         tile_size = size;
 472     else
 473         tile_num = ((size.width + tile_size.width - 1)/tile_size.width)*
 474                    ((size.height + tile_size.height - 1)/tile_size.height);
 475
 476     CV_CALL( tiles = (double*)cvAlloc( tile_num*10*sizeof(double)));
 477
 478     for( y = 0, k = 0; y < size.height; y += tile_size.height )
 479     {
 480         CvSize cur_tile_size = tile_size;
 481         if( y + cur_tile_size.height > size.height )
 482             cur_tile_size.height = size.height - y;
 483
 484         for( x = 0; x < size.width; x += tile_size.width, k++ )
 485         {
 486             if( x + cur_tile_size.width > size.width )
 487                 cur_tile_size.width = size.width - x;
 488
 489             assert( k < tile_num );
 490
 491             IPPI_CALL( func( mat->data.ptr + y*mat->step + x*pix_size,
 492                              mat->step, cur_tile_size, cn, coi, tiles + k*10 ));
 493         }
 494     }
 495
 496     icvAccumulateMoments( tiles, size, tile_size, moments );
 497
 498     __END__;
 499
 500     if( ippmomentstate )
 501         icvMomentFree_64f_p( ippmomentstate );
 502
 503     cvFree( &tiles );
 504 }
 505
 506 /*F///////////////////////////////////////////////////////////////////////////////////////
 507 //    Name: cvGetHuMoments
 508 //    Purpose: Returns Hu moments
 509 //    Context:
 510 //    Parameters:
 511 //      mState  - moment structure filled by one of the icvMoments[Binary]*** function
 512 //      HuState - pointer to output structure containing seven Hu moments
 513 //    Returns:
 514 //      CV_NO_ERR if success or error code
 515 //    Notes:
 516 //F*/
 517 CV_IMPL void
 518 cvGetHuMoments( CvMoments * mState, CvHuMoments * HuState )
 519 {
 520     CV_FUNCNAME( "cvGetHuMoments" );
 521
 522     __BEGIN__;
 523
 524     if( !mState || !HuState )
 525         CV_ERROR_FROM_STATUS( CV_NULLPTR_ERR );
 526
 527     {
 528         double m00s = mState->inv_sqrt_m00, m00 = m00s * m00s, s2 = m00 * m00, s3 = s2 * m00s;
 529
 530         double nu20 = mState->mu20 * s2,
 531             nu11 = mState->mu11 * s2,
 532             nu02 = mState->mu02 * s2,
 533             nu30 = mState->mu30 * s3,
 534             nu21 = mState->mu21 * s3, nu12 = mState->mu12 * s3, nu03 = mState->mu03 * s3;
 535
 536         double t0 = nu30 + nu12;
 537         double t1 = nu21 + nu03;
 538
 539         double q0 = t0 * t0, q1 = t1 * t1;
 540
 541         double n4 = 4 * nu11;
 542         double s = nu20 + nu02;
 543         double d = nu20 - nu02;
 544
 545         HuState->hu1 = s;
 546         HuState->hu2 = d * d + n4 * nu11;
 547         HuState->hu4 = q0 + q1;
 548         HuState->hu6 = d * (q0 - q1) + n4 * t0 * t1;
 549
 550         t0 *= q0 - 3 * q1;
 551         t1 *= 3 * q0 - q1;
 552
 553         q0 = nu30 - 3 * nu12;
 554         q1 = 3 * nu21 - nu03;
 555
 556         HuState->hu3 = q0 * q0 + q1 * q1;
 557         HuState->hu5 = q0 * t0 + q1 * t1;
 558         HuState->hu7 = q1 * t0 - q0 * t1;
 559     }
 560
 561     __END__;
 562 }
 563
 564
 565 /*F///////////////////////////////////////////////////////////////////////////////////////
 566 //    Name: cvGetSpatialMoment
 567 //    Purpose:  Returns spatial moment(x_order, y_order) which is determined as:
 568 //              m(x_o,y_o) = sum (x ^ x_o)*(y ^ y_o)*I(x,y)
 569 //              0 <= x_o, y_o; x_o + y_o <= 3
 570 //    Context:
 571 //    Parameters:
 572 //      mom  - moment structure filled by one of the icvMoments[Binary]*** function
 573 //      x_order - x order of the moment
 574 //      y_order - y order of the moment
 575 //    Returns:
 576 //      moment value or large negative number (-DBL_MAX) if error
 577 //    Notes:
 578 //F*/
 579 CV_IMPL double
 580 cvGetSpatialMoment( CvMoments * moments, int x_order, int y_order )
 581 {
 582     int order = x_order + y_order;
 583     double moment = -DBL_MAX;
 584
 585     CV_FUNCNAME( "cvGetSpatialMoment" );
 586
 587     __BEGIN__;
 588
 589     if( !moments )
 590         CV_ERROR_FROM_STATUS( CV_NULLPTR_ERR );
 591     if( (x_order | y_order) < 0 || order > 3 )
 592         CV_ERROR_FROM_STATUS( CV_BADRANGE_ERR );
 593
 594     moment = (&(moments->m00))[order + (order >> 1) + (order > 2) * 2 + y_order];
 595
 596     __END__;
 597
 598     return moment;
 599 }
 600
 601
 602 /*F///////////////////////////////////////////////////////////////////////////////////////
 603 //    Name: cvGetCentralMoment
 604 //    Purpose:  Returns central moment(x_order, y_order) which is determined as:
 605 //              mu(x_o,y_o) = sum ((x - xc)^ x_o)*((y - yc) ^ y_o)*I(x,y)
 606 //              0 <= x_o, y_o; x_o + y_o <= 3,
 607 //              (xc, yc) = (m10/m00,m01/m00) - center of gravity
 608 //    Context:
 609 //    Parameters:
 610 //      mom  - moment structure filled by one of the icvMoments[Binary]*** function
 611 //      x_order - x order of the moment
 612 //      y_order - y order of the moment
 613 //    Returns:
 614 //      moment value or large negative number (-DBL_MAX) if error
 615 //    Notes:
 616 //F*/
 617 CV_IMPL double
 618 cvGetCentralMoment( CvMoments * moments, int x_order, int y_order )
 619 {
 620     int order = x_order + y_order;
 621     double mu = 0;
 622
 623     CV_FUNCNAME( "cvGetCentralMoment" );
 624
 625     __BEGIN__;
 626
 627     if( !moments )
 628         CV_ERROR_FROM_STATUS( CV_NULLPTR_ERR );
 629     if( (x_order | y_order) < 0 || order > 3 )
 630         CV_ERROR_FROM_STATUS( CV_BADRANGE_ERR );
 631
 632     if( order >= 2 )
 633     {
 634         mu = (&(moments->m00))[4 + order * 3 + y_order];
 635     }
 636     else if( order == 0 )
 637         mu = moments->m00;
 638
 639     __END__;
 640
 641     return mu;
 642 }
 643
 644
 645 /*F///////////////////////////////////////////////////////////////////////////////////////
 646 //    Name: cvGetNormalizedCentralMoment
 647 //    Purpose: Returns normalized central moment(x_order,y_order) which is determined as:
 648 //             nu(x_o,y_o) = mu(x_o, y_o)/(m00 ^ (((x_o + y_o)/2) + 1))
 649 //             0 <= x_o, y_o; x_o + y_o <= 3,
 650 //             (xc, yc) = (m10/m00,m01/m00) - center of gravity
 651 //    Context:
 652 //    Parameters:
 653 //      mom  - moment structure filled by one of the icvMoments[Binary]*** function
 654 //      x_order - x order of the moment
 655 //      y_order - y order of the moment
 656 //    Returns:
 657 //      moment value or large negative number (-DBL_MAX) if error
 658 //    Notes:
 659 //F*/
 660 CV_IMPL double
 661 cvGetNormalizedCentralMoment( CvMoments * moments, int x_order, int y_order )
 662 {
 663     int order = x_order + y_order;
 664     double mu = 0;
 665     double m00s, m00;
 666
 667     CV_FUNCNAME( "cvGetCentralNormalizedMoment" );
 668
 669     __BEGIN__;
 670
 671     mu = cvGetCentralMoment( moments, x_order, y_order );
 672     CV_CHECK();
 673
 674     m00s = moments->inv_sqrt_m00;
 675     m00 = m00s * m00s;
 676
 677     while( --order >= 0 )
 678         m00 *= m00s;
 679     mu *= m00;
 680
 681     __END__;
 682
 683     return mu;
 684 }
 685
 686
 687 /* End of file. */