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