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 icvCompleteMomentState(CvMoments * moments)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 icvContourMoments(CvSeq * contour,CvMoments * moments)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 icvAccumulateMoments(double * tiles,CvSize size,CvSize tile_size,CvMoments * moments)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 cvMoments(const void * array,CvMoments * moments,int binary)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 cvGetHuMoments(CvMoments * mState,CvHuMoments * HuState)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 cvGetSpatialMoment(CvMoments * moments,int x_order,int y_order)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 cvGetCentralMoment(CvMoments * moments,int x_order,int y_order)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 cvGetNormalizedCentralMoment(CvMoments * moments,int x_order,int y_order)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. */ 688